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INTRODUCTION 


The  family  of  computer  programs  described  in  this  report  all  calcu- 
late a six  degree  of  freedom  simulation  of  the  flight  of  a rigid  body  through 
the  earth's  atmosphere.  The  elements  of  the  transformation  matrix  from 
inertial  axis  to  body  axis  are  produced  by  direct  integration  rather  than 
through  Euler  angles.  Aerodynamic  force  and  moment  coefficients  are  in- 
troduced as  tables  as  a function  of  Mach  number,  total  angle  of  attack,  and 
roll  position  of  the  plane  of  angle  of  attack.  They  are  so  written  as  to  make 
the  use  of  body-fixed  axes  possible  even  when  the  integration  proceeds  in 
non-rolling  coordinates.  The  vehicle  is  allowed  an  initial  coast,  a powered 
phase,  and  a final  coasting  phase. 

The  family  of  programs  includes: 

PANOL  2 - for  symmetric  missiles  and  projectiles  where 
products  of  inertia  vanish  and  the  center  of  mass  is  on  the  longitudinal  axis. 

PANOL  3 - for  unsymmetric  missiles  and  projectiles  where 
products  of  inertia  may  appear  and  the  center  of  mass  may  lie  off  the  longi- 
tudinal axis,  but  the  mass  is  constant. 

PANOL  4 - for  yaw  sonde  simulations  of  symmetric  projectiles . 

PANOL  5 - for  guided  missiles  and  projectiles  with  a target 
consisting  of  an  upright,  rectangular  plane  subdivided  into  rectangles  of 
varying  brightness;  with  thrust  control;  on  a non-rotating  spherical  earth. 


The  computer  program  input  instructions  are  included  as  Appendix  A . 
Output  formats,  titles,  and  units  are  contained  in  Appendix  B. 

MATHEMATICAL  FORMULATION  OF  PROBLEM 


Rigid  Body  Equations  of  Motion 

The  force  and  moment  equations  can  be  written  in  the  form 

(1) 
(2) 


F = / a dm 
M = / r X a dm 


J 

1 


) 

\ 


s 
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These  integrals  extend  over  all  mass  particles  within  the  geometric  bound- 
ary of  the  vehicle  at  a given  instant  of  time.  This  boundary  lies  on  the 
vehicle  surface  and  includes  the  rocket  nozzle  exit  surface. 

The  absolute  acceleration  of  a mass  particle  within  the  vehicle's 
geometric  boundary  is 

a = a +05xr  + o>x(coxr)+2(joxV  +a  (3) 

o r r 


The  mass  contained  within  the  geometric  boundary  of  the  vehicle 
at  any  instant  is  composed  of  the  mass  of  the  vehicle  structure  and  the 
mass  of  the  unburned  propellant,  the  sum  being  designated  as  m^,  plus 
the  mass  (nig)  of  the  propellant  gases  which  have  not  yet  crossed  the 
nozzle  exit  plane.  Accordingly,  Equations  1 and  2 can  be  rewritten  as 


F = / a dm  + / a dm 
m m 

s g 


M = / r X a dm  + / r x a dm 


m 


ro 


g 


(4) 

(5) 


For  the  time  being  the  computational  axes  will  be  rigidly  fixed  to 
the  ■ ■ ’cle  so  that  in  the  first  integral  appearing  in  Equations  4 and  5, 


V 


o.  Also,  since  m »m  , the  integrals 
S fi 


/ [a  +a)xr  + a3x(o3xr)]dm 
m o 
g 


and 


/ rx[a  +a3xr  + mx(a5xr)]  dm 
m o 

g 


will  be  neglected.  Subject  to  these  restrictions.  Equations  4 and  5 can  be 
written  as 


F = / [a  +o3xr  + cox(coxr)]  dm  + 2/  oi  x V dm  + 
mo  m r 

s g 


/ a dm 
m r 
g 


(6) 


2 


I 


and 


M=/  rxia  +aixr  + o3x((joxr)l  dm  + 
m o 

s 


2/  rx(coxV)dm  + / rxa  dm 
m r m r 

g g 


(7) 


The  first  of  the  three  integrals  in  Equations  6 and  7 will  be  referred 
to  as  rigid  body  integrals;  the  second,  as  jet  damping  integrals;  and  the 
third,  as  thrust  integrals. 


Rigid  Body  Integrals 


Noting  that  a^,  to,  and  co  are  common  to  all  mass  particles 


comprising  m^,  the  first  integral  in  Equation  6 can  be  written  as 


A =a  m +cox/  rdm+mlmx/^  r dm) 
r o s m m 

s s 


(8) 


The  integral  in  Equation  8 is  equal  to  m^  r^  where  the 

vector  r describes  the  location  of  the  center  of  gravity  of  the  vehicle 
c 


with  respect  to  the  origin  of  the  body  axes.  The  center  of  gravity  of  the 
vehicle  will  be  restricted  to  lie  on  the  axis  so  the  r^  = x^i . Equation  8 
can  be  expanded  to  read 

= m^  [ ^ - Cq^  + r^ ) ] i + m^  [ a^  + (r  + p q)  X^]  j + 

(9) 


m [a  - (q  - pr)  X ] k 
so  c 

z 


The  first  integral  in  Equation  7 can  be  written  as 


G =m  r xa  + / rx(o5xr)dm  + / rx[a5x(o3xr)]dm 
r s c o m ni 

s s 


(10) 


Imposing  the  requirement  that  the  body  axes  are  principal  axes  of  inertia. 
Equation  10  can  be  expanded  to  yield 


G = [p  I + rq  (I  - I ) 1 
r X ^ z y 

s s s 


[ql  -rp  (I^  - ) -m^  ] ; 


[rl  + pq  (I  ~ I ) + m X a ]k 
•>  r'l',,  V SCO 


y X 

s s 


(11) 


In  some  versions  of  this  program  I = I.  : however,  since 

^s  ^s 

they  are  sometimes  distinct,  they  will  be  designated  separately.  In  some 
versions  of  this  program  the  body  axes  are  not  principal  axes,  but  in  those 
cases,  thrust  must  be  zero  and  mass  is  constant. 


Jet  Damping  Integrals 

These  integrals  are  the  negative  of  the  force  and  moment  re- 
actions applied  to  the  vehicle  by  the  Coriolis  acceleration  of  the  propellant 
gases  prior  to  passage  through  the  nozzle  exit  plane. 

The  jet  damping  force  is  related  to  the  second  integral  in 
Equation  6 by  the  expression 

= -2/  (oxVdm  = 2[/  (rV  -qV  ) dm]  i + 

J m r m r ^ r 

g g y z 

2 [/  (pV  -rV  )dm]t+2[/  (q  - p ) dm]  k 
g z ^x  g X y 

(12) 


If  the  velocity  field  of  the  propellant  gases  is  axisymmetric  with  respect 
to  the  axis,  then  terms  involving  and  drop  out  by  cancellation 

^y  '^z 

when  we  integrate  over  opposite  sides  of  the  axis  and 


where 


Fj  = - J,  (rj  -qk) 


J.  = 2/  V dm 
' m r 
g X 


(13) 
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Similarly,  the  jet  damping  moment  is  given  by  the  expression 


M^  = -2/  rx(o)xV)dm 
J m r 

g 

= 2[f  (qyV  +rzV  -pyV  -pzV  ) dm]  i + 
m r r r r 

g X X y z 


2[/  (rzV  +pxV  -qzV  -qxV  )dm]  j + 
g y y z X 

2[f  (pxV  +qyV  -rxV  - ryV  )dm]  k (14) 
m *^  r ^■'r  r •'r 

g z z X y 

In  addition  to  the  assumption  that  the  propellant  gas  flow  is  symmetric 
about  the  axis,  it  is  also  assumed  that  the  passage  through  which  these 

gases  flow  has  a sufficiently  small  diameter  and  curvature  such  that  the 
approximations  x»y , x»z , V »V  , and  V »V  can  be  made  in 
Equation  14.  ^y 


Therefore, 


where 


Mj  = - Jj  (qj  + r k) 


J,  = 2/  X V dm 


Thrust  Integrals 

The  thrust  integrals  are  the  negative  of  the  thrust  forces  and 
moments  produced  by  the  motor . These  forces  and  moments  account  for 
the  reactions  due  to  the  momentum  flux  of  the  propellant  gases  passing 
through  the  nozzle  exit  plane,  plus  the  pressure  acting  on  the  nozzle  exit 


The  thrust  force  is 


F = / a dm 
T m r 


In  general 


1 


a(p. 


r a dm  = / 

m r m p 

g g g 


V ) 
r 


at 


dm  + / 


V - 
e 


n )da 


where  the  first  integral  on  the  right  side  of  this  expression  is  a volume 
integral  similar  to  others  previously  treated.  The  second  integral  is  a 
surface  integral  evaluated  over  the  nozzle  exit  surface.  If  the  propellant 
gas  flow  is  steady  then  the  thrust  force  can  be  shown  to  be 

e 


where 


and 


n = (p  - p)  n 
e e e 


p = -p  { V • n ) 
e e e e 


p , p , and  V are  the  static  pressure,  density  and  velocity  of  the  pro- 
6.  6 6 

pellant  gases  at  the  nozzle  exit  plane:  n^  is  the  unit  inward  normal  to  this 

plane;  and  p is  the  ambient  atmospheric  pressure.  Assuming  that  p , p , 
-1-  -V  5.  ® 

V , and  n are  constant  over  the  nozzle  exit  plane,  and  that  V and  n are 
0 0 ^ 0 0 
parallel  at  all  points  on  this  plane,  then  Equation  17  becomes 


F + (p  - p)  A n + V / p (V  . n )da 
T *^e  ^ e e e A ^e  e e 

e 


Now  note  that 


/.  p (V  »n)da  = p (V  ‘n)/.  da  = m<o 
A e e e e e e A 

e e 


and  also  that 


V = - V n 
e e e 


so  that 


Tn 


e 


6 


where 


T = [(Pg  - p)A^  - mVg] 

The  orientation  of  the  unit  vector  n with  respect  to  the  body  frame  is 
shown  below: 

I ® 


The  vector  n^  in  the  body  frame  can  be  written  as 

n = cos (62 ) cos (Ej ) i + cos  (£2 ) sin (Ej ) j + sin  (E2 ) k 
3 

If  the  angles  Ej  and  E2  are  sufficiently  small  then 

= T ( i + El  j + E2  k)  (18) 

The  thrust  moment  is 

M = -/  r X a dm  (19) 

T m r 

g 


7 


I 


r 


m 


For  the  case  where  the  propellant  gas  flow  is  steady,  similar  to  the  case  of 
thrust  force,  Equation  19  can  be  written  as 


(20) 


e e e 


where  r^  is  the  position  vector  of  any  point  on  the  exhaust  exit  plane. 
Under  the  assumptions  previously  stated,  the  thrust  moment  can  be  written 
as 

M = [p  - p - p V ( V • n ) ] / ( r X n ) da 

T e ee  e e A e e 

e 


or  simply 


M = (p  - p + p )/  (r  X n ) da 
-r  e A e e 

e 


(21) 


The  integral  in  Equation  21  is  readily  expanded  to  yield 


U [(ye  £2  - Zg  El)  i + 

e e 


(^e  ■ ^ [^e  " ^e^  ^ ^ 


resulting  in  the  expression 


M.J,  = T [ (sj  y^  - Si  2g)  i + (Zg  - x^  Ej ) j + ( Ej  - yg)  k ] (22) 


where 

/.  X da  = X A , etc 
A e e e 
e 

so  that  X , y , 2 are  the  coordinates  of  the  centroid  of  A measured  in 
e -^e  e e 

the  body  frame. 
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Gravity  Forces  and  Moments 


The  gravity  force  is  given  by  the  expression 


where  Ris  the  position  vector  of  an  infinitesmal  mass  (dm)  of  the  vehicle 
with  respect  to  the  Earth's  center.  For  a vehicle  of  reasonably  small  di- 
mensions, flying  in  the  atmosphere,  a good  approximation  is 

R * R 

where  R is  the  position  vector  of  the  center  of  gravity  of  the  vehicle  with 
respect  to  the  Earth's  center.  Neglecting  the  effect  of  gravity  on  the  pro- 
pellant gas  particles,  then 

F = — £ (X  t+y  t+z  k)  (24) 

G ,-,3  o ^o^  o 


where  x^,  y^,  are  the  components  of  R in  the  body  frame. 

The  gravity  moment  is  given  by  the  expression 


M„=  gR'/^dm 

G o e I B I , 


Again,  assuming  that  R * R Equation  25  becomes 


®o  ^E  ^ ^ 

= Rx/rdm 


Neglecting  the  mass  of  the  propellant  gases  and  recalling  that 


rdm  = m r =m  x i 
s c s c 


the  gravity  moment  can  be  written  as 


■"s  8„ 


Aerodynamic  Forces  and  Moments 

The  aerodynamic  forces  and  moments  acting  on  the  vehicle 
can  be  written  as 


= QA  i ♦ Cy  J * Cj  k ) 


Sf^  = QAd  (Cj_r*C|^t+Cj,k)  (28) 

These  equations  will  be  developed  in  more  detail  in  a subsequent  section. 

—Ik 

Substituting  the  definitions  of  A from  Equation  8,  F from 

Equation  12.  F.^  from  Equation  16.  F^  from  Equation  23  and  F^  from  Equa- 
tion 27  into  Equation  6 we  get 

A = F ^ + F ^ + F + F y 
r G A T J 

Similarly  the  definitions  of  G from  Equation  10,  from 
— i —X  r J 

Equation  14.  M from  Equation  19,  M_  from  Equation  25  and  M.  from 
1 G A 

Equation  28 , substituted  into  Equation  7 yield 

G = + M . + N?  + M y 

r G A T J 

Expanding  the  forces  and  accelerations  above  and  rearranging 
the  terms , we  get 


' . . I,  . i.*  . , ... 


The  force  Equation  29  will  now  be  written  in  terms  of  inertial 
frame  components . The  transformation  matrix  [£  ] transforms  vector  com- 

ponents  in  the  body  frame  to  vector  components  in  the  inertial  frame.  The 
components  of  the  vectors  a^  and  R in  these  two  frames  are,  therefore,  re- 
lated by  the  expressions 


a 

o 

X 

1 

1 

a 

U 

y 

a 

• • 

0 

z 

^R 

and 


X 

0 

R 

^O 

^R 

z 

0 

R 

Note  that  a^  was  the  acceleration  of  the  origin  of  the  body  frame 

with  respect  to  an  inertial  coordinate  system  instantaneously  parallel  with 
the  body  frame. 


Substituting  these  expressions  into  Equation  29  and  premulti- 
plying each  side  by  [£„]  , the  force  equation  becomes 


1 — 

IX 

•• 

g Ri, 
“0  E 

^R 

1 rT 

^R 

• • 

Z 

R 

f 

F 

X 

1 

’o  1 

q*  + r^ 

F 

y 

T 

+ — 
m 

s 

El 

Ji 

"m 

s 

r 

+ x 

c 

-i  - pq 

1 

F 

z 

-q 

q - pr 

L J 

_ j 

- 

> 
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Equations  30  and  31  now  will  be  written  in  a form  which  will  allow  the  use 
of  a body  frame  rigidly  attached  to  the  vehicle,  for  which  they  are  already 
valid,  or  for  a body  frame  which  does  not  roll  with  the  missile.  The  neces- 
sary transformations  of  the  Force  Equation  31  and  the  Moment  Equation  30 
will  be  presented  next. 

Symbolically.  Equation  31  can  be  written  as 

^ ^ ^ 

p = F^  + fC  ]{F  +F  +F  +F) 

G ^ ' A T J r 

The  vectors  P and  F^,  which  describe  the  inertial  acceleration  of  the 
body  axes  origin  and  the  gravity  forces,  are  already  in  the  desired  form 
since  they  are  described  in  terms  of  variables  resolved  in  the  inertial 
frame . 

The  matrix  [£g]  can  be  calculated  so  that  it  will  transform 
a vector  in  either  a rolling  or  nonrolling  body  frame  into  a vector  in  the 
inertial  frame.  This  is  shown  in  Reference  1 and  repeated  here  for  con- 
venience. The  [2g]  matrix  is  calculated  by  integrating  the  set  of  nine 
differential  equations 


[Kg]  = [Kg]  [n] 


where 


m = 


-q  kp  o 


If  k = 1 the  body  axes  are  rolling,  and  if  k = 0 they  are  nonrolling.  A cor- 
rection procedure  for  the  minimization  of  the  nonorthogonality  introduced 
by  the  numerical  integration  of  these  equations  is  incorporated . 


The  calculation  of  the  aerodynamic  force  vector  F , as  will 

A A 

be  seen  later,  is  accomplished  in  such  a way  that  either  the  rolling  or 
nonrolling  computational  frames  can  be  used.  The  thrust  jet  damping 
Fj  and  inertia  Fj  force  vectors  now  must  be  transformed. 


I 


The  orientation  of  the  rolling  frame  (x,  y,  z)  with  respect 
to  the  rolling  or  nonrolling  computational  frame  (x',  y',  x')  is  depicted 
in  the  sketch  below. 


The  components  of  any  vector  in  the  rolling  or  nonrolling  computational 
frame  are  related  to  its  components  in  the  rolling . body-fixed  frame  by 
the  expression 

V = [£(p']  V ’ 

where  primes  denote  vectors  resolved  in  the  computational  frame. 


14 


The  components  of  the  time  rate  of  change  of  V are  given  by 
Ik  * _i.  » Ij.  » 

V = [e<p']  V + [Cep']  V 


where 


[Cep']  = 


o cos  ep'  sin  ep' 

o -sin  ep'  cos  ep' 


Using  these  relations,  the  vectors  F , F anel  F transform 

1 u 1 


» T T 

Ft  Ft  = - 


^ r T 1 

F^  = [Ceo]^  = m 

J J m 


Fj'  = [2^>']^F^= 


El  cos(ep')  - E2  sin(ep') 

£2  cos(ep')  + El  sin(q)')  j 
o I v;here 

r'  ; , r'  = r cos(ep')  + q sin(ep') 

q'  = q cos(ep')  - r sin(ep') 


anel 


q 

q' 


>2  + J.I2  ~j 


where 


anel 


equation 


r'  = [r  + q (l-k)p]cos  (cp')  + [q-r  (l-k)p]  sin  (to') 

q' = [q-r  (l-k)p]cos  (ep')  - [r+q  (l-k)p]  sin  (ep') 

Since  the  roll  angle  ep'  is  determined  from  the  differential 


with  initial  condition 


ep’  = (l-k)p 
ep'  (o)  = o 
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T" 


Equation  31  can  now  be  written  in  a form  which  applies  to  either  a rolling 
or  nonrolling  computational  frame,  depending  upon  the  choice  of  k.  That  is, 


0>  «• 

- 

( 

-| 

^R 

F i 
X 1 

Y 

g R* 
o E 

1 

I 

F ! 

R 

1 R 1 ’ 

R 

B ' 

m 

s 

y ' 

F 1 

R 

R 

z 

_ _ 

^ 

_ J 

1 

o 

1 

>.r^  - 

\ 

F-i  cos((p')  - S2  sin(cp') 

J, 

m 

r 

+ X 

c 

-f  - kpq 

s 

i 

Cj  cos(cp')  + El  sin(cp') 

-q 

q - kpr 

, 1 

(36) 


The  primes  on  the  angular  velocities  are  dropped  since  it  is  understood 
that  these  variables  are  measured  in  the  particular  frame  which  is  being 
used  for  computation. 

The  transformation  of  Equation  30  is  accomplished  in  a 
similar  fashion;  the  symbolic  form  of  this  equation  is 

H = M ^ + M ^ + Mj  + 


Writing  the  rate  of  change  of  angular  momentum  in  the 
computational  frame,  which  may  be  rolling  or  nonrolling,  we  get 


" — 

P I 

X 

P ( 

X 

S 

s 

+ [O] 

q I , 

^s 

y 

s 

r I 

r I 

z 

z 

s 

L ^ . 
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It  should  be  pointed  out  that  this  form  is  good  only  for  those  cases  where 


1 = I . The  vector  H could  easily  be  put  in  a form  compatible  with 


both  the  rolling  and  nonrolling  frames  for  the  case  where  I 1 , but 

^s  s 


some  of  the  advantage  of  using  the  nonrolling  frame  would  be  lost,  since 
the  moment  of  inertia  matrix  would  not  be  constant  in  the  nonrolling  com- 
putational frame  and  an  additional  term  would  have  to  be  calculated.  Fur- 
thermore , the  computation  of  the  moment  of  inertia  and  its  rate  of  change 
would  then  require  an  accurate  knowledge  of  the  roll  position  of  the  body 
fixed  axes,  precisely  what  we  are  trying  to  avoid  computing  by  intro- 
ducing the  nonrolling  computational  axes. 


As  in  the  case  of  the  aerodynamic  forces,  the  vector  will 


be  calculated  in  a manner  compatible  with  both  the  rolling  and  nonrolling 
frames.  The  gravity,  thrust,  jet  damping, and  inertia  moment  vectors  are 
transformed  as  were  the  force  vectors. 


Using  Equation  33,  the  transformed  vectors  are 


M, 


mg  X 

so  E c 


R I 3 


~y' 


I 

=T 


£2 


e,  7 


..e 


(z  - Ej  X ) cos(cp')  - (Ej  X - y ) sin(cp') 


(ej  X - y ) cos(tp')  + (z  - E,  X ) sin{cp'] 
e e 0*0  T 


M. 


- J, 


o 

q' 

r' 


1 


M,  = m X 

I s c 


yj 
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( 


so  that  the  r omeiit  equation  becomes 


r 

— 

l‘x 

S 

pi 

^ X 

s 

0 

i 

S’ 

s 

= - [n] 

qi 

y 

s 

, "■s8o'’e  '■0 
^ 1 R|’ 

z 

0 

f 1 

L 

rl 

z 

s 

0 

1 

I 

M 


M 


M 


+ T 


E.  Yg  - El  Z3 


(z  - £2  X ) cos(cp')  - (El  X - y ) sin(cp') 
^ 6 6 


(Ej  X - y ) cos((p')  + (z  - £2  X ) sin(cp') 


e e 


- J, 


— — 

— 

0 

0 

q 

+ m X 

a 

s c 

0 

z 

r 

-a 

0 

y 



(37) 


Differential  Equations  36  and  37  describe  the  motion  of  a 

rolling  or  nonrolling  computational  frame  which  has  its  origin  fixed  at 

the  burnout  center  of  gravity  of  the  vehicle.  These  equations  are  not  in 

a form  suitable  for  integration  because  of  the  terms  f and  q in  Equation  36, 

and  a and  a in  Equation  37. 
o o 

y X 

At  this  point  it  should  be  recalled  that  in  order  to  calculate 
the  aerodynamic  forces  and  moments,  the  velocity  of  the  vehicle  structure 
with  respect  to  the  local  air  mass  must  be  known.  However,  the  velocity 
of  the  origin  of  the  body  fram.e  is  not  the  same  as  this  velocity  because  of 
the  rotational  acceleration  components  which  account  for  the  fact  that  the 
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burnout  center  of  gravity  position  does  not  coincide  with  the  instantaneous 
center  of  gravity  position.  If  the  last  term  in  Equation  36  is  suppressed, 
then  these  equations  describe  the  position  of  the  instantaneous  center  of 
gravity  and  the  velocity  of  a fixed  point  coinciding  with  the  instantaneous 
center  of  gravity  at  any  given  time. 

The  Moment  Equation,  37,  can  be  readily  put  into  a form  suit- 
able for  integration  if  Equation  29  is  used  to  eliminate  a and  a^  and  if 

y 2 

the  remaining  inertia  terms  are  transformed  so  that  they  are  compatible 
with  both  computational  frames. 


final  form: 


The  force  and  moment  equations  can  now  be  written,  in  their 


R 


1 nl 

L .1 


f”  “I 

' X 
I R. 


= - ®o'>E 


M 


R 


i ’^Ri  ^~B 


R i 

J 


r 

j 

1 

1 " 
! ° 

j 

' El  COS(cp') 

- Ej  Sin  (<p') 

m ' 

Ej  cos  (9') 

+ Ej  sin  (cp') 

S ' 

-q 

L- 

_i 

_ J 

(38) 
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I 


'-2 


e,  z 


T 1 (z  - Eo  X ) cos((p')  - (e,  X - y ) sin((p') 

10"*  P ^ ‘p'o  ^ 


e ' e 


i (e,  X - y J cos  (<f ')  + (z  - Ez  X ) sin  Icp') 

] G Q 0 0 


- J, 


(39) 


Since  Equation  38  can  be  simply  transformed  to  a set  of  six  first 
order  differential  equations.  Equations  32,  35,  38,  and  39  are  equivalent 
to  a set  of  nineteen  first-order  differential  equations  describing  the  ve- 
hicle's motion . 

Equations  Describing  Physical  Characteristics  of  Vehicle 

The  physical  characteristics  of  the  vehicle  are  calculated  based  on 
two  models  of  the  rocket  motor.  Both  have  a single  period  during  which  thrust 
occurs,  and  both  have  constant  mass  flow  rate  during  burning  with  igni- 
tion occurring  at  t = tj  and  burnout  occurring  at  t = t2  . Two  sets  of  equa- 
tions are  presented;  one  for  an  end  burning  grain,  and  the  other  for  a 
radially  burning  grain. 

The  calculations  of  thrust  and  mass  are  identical  for  both  cases.  The 
thrust  equation  is 


T = [(p  - p)  A 
e e 


-mV] 

e 


All  quantities  except  the  ambient  atmospheric  pressure  (p)  are  input  con- 
stants. The  ambient  pressure  is  a function  of  altitude  as  given  in  Reference  3. 
The  instantaneous  vehicle  mass  and  propellant  mass  are 


and 


m = m + m 
s p 


+ m (t  - tj ) 
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All  quantities  except  t are  input  constants,  with  m beinq  a negative  quan- 
tity . 


The  remainder  of  the  calculations  depend  on  the  mode  of  burning 
and  will  be  considered  in  two  parts.  Refer  to  Figure  1 for  a sketch  showing 
the  dimensions  involved. 

End  Burning  Crain 

1.  Center  of  gravity  location: 


m 


X = 

c m 


2.  Axial  moment  of  inertia: 


2 


m r 

I = I + 

X X 2 
s 


Transverse  moment  of  inertia: 


I = I = I + m 

ys  ^ y p 


Li!  ^ - Xg  ) ^ ^ ^ 2 


12 


+ X, 


4.  Jet  damping  factors: 


2/  V dm  = 2// 
m r 


n V da  dx 

^ X 


■ 

Ji  = 2mj  dx  = 2m(x3  - X4 ) 


Similarly 

J2  = 2/  X V dm  = m (X3*  - X4*) 
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For  this  case.  Xj  , , and  rj  are  input  constants  and 


X2  = r (x,  + x,: 


Xj  = Xj 


m (t  - t, ) 


O r»  Pr 


Radially  Burning  Crain 

1.  Center  of  gravity  location: 


m 

X = X2 

C 


2.  Axial  moment  of  inertia: 


I = I + ^ (n*  +r2*) 

X X 2 ' ‘ 

s 


Transverse  moment  of  inertia: 


Zs  y p 4 


4.  Jet  damping  factors: 


Ji  = 2//p  V da  dx 
R r 
® X 


In  this  case,  however, 


/p  V da 
r 


Xi  > X > X3 


X,  > X > Xi 


/ • ».  >.1.  - . - . 


So  that 


X,  Xj 

J,  = 2/  m ( ^ ) dx  + 2/  m dx 

X,  ^3-X, 


J,  = m (Xi  + Xj  - 2 X4 ) 


Similarly , 


J2=2//xp  V dadx 

E r 
® X 

J2  = ^3^  + Xi  Xj  - 3 X4*) 


For  this  case  Xx  , X2  = X2  , X3  = X3  , X4  , and  rj  are  input 
constants  and  ° o 


1^2  = 


■:  I 

m (t  - ti ) 

Pp  (Xi  - X3) 


Reference  5 contains  additional  information  on  the  calculation  of  the  jet 
damping  factors. 


AERODYNAMIC  FORCE  AND  MOMENT  EQUATIONS 


The  aerodynamic  forces  and  moments  as  input  to  the  program  are 
always  calculated  as  originating  in  their  proper  coordinate  system;  i.e. 


the  trim  normal  force,  F 


is  resolved  into  body  fixed  coordinates  no 


matter  whether  integration  takes  place  in  body  fixed  or  nonrolling  coor- 
dinates. The  computational  frame,  subscripted  AB,  into  which  all  aero- 
dynamic forces  and  moments  are  resolved  is  shown  on  Figure  2.  The 
following  roll  position  and  roll  relationship  angles  are  needed  for  the 
analysis: 


0 Azimuth  angle  of  the  angle  of  attack  plane  in  the  AB  frame 
qi'  Roll  angle  between  AB  frame  and  rolling,  body-fixed  frame 
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I 


J 


Azimuth  angle  of  the  angle  of  attack  plane  in  the  rolling, 
body-fixed  plane 


(i ' Reduced  <p;  it  is  [<p  Modulo  !^]  where  % is  the  repeat 
roll  angle  of  the  aerodynamic  symmetry 

The  total  aerodynamic  force  and  moment  coefficients,  in  terms  of 
the  Taylor  series  expansion  coefficients,  are 


C„  = C + C 
X X X 2V . 

P A 


(for  a normal  missile  < o) 


C =C  cos(«i)  +C  sin(0)  +C  cos(ci)  + 

Y - y z yp  2V^ 

C [6„  cos(<p')  - 6.  sin(cp')] 

z_  B ^ A 

^A 

(for  a normal  missile 
if  (6  < 90°;  Cy  < o) 

C_  = C cos(d)  - C sin(d)  - C,,  (?n — ) sin  (ti)  + 

Z z y Yp  2V^ 

C [6.  cos((fl')  +6  sin(cp')] 

z_  A B 

®A 

(for  a normal  missile 

= o by  the  axisymetric  assumption) 


Cl.c,*c  a*c,  (f|-) 

6 P A 


C,,  E C cos(«i)  + C sin({S)  + C (^^ — ) sin  («i)  + 

Mm  n n^  2V^ 

q A 6 


(for  a statically  stable  missile. 


C = < o) 
m 


( 


T 


C = C cos  (0  -C  sin  (0  +C  — cos(0)  *C  — 

\ n m n 2\\  m 2'  , 

p A q A 


[6gCos((p')  6^  sin  (cf/)  ] (■}( 

(for  a normal  missile, 

= 0 by  the  axisymm.etric 
assumption! 


where  C , C 

X X 


1 , C , C , C , C , and  C 

Zr  I m n m. 

6 6 p q p 6, 


are  functions 


p p t)^  6 p q p 

of  Mach  number  and  total  angle  of  attack.  The  coefficients  C . C , C , 

y z i 

and  are  functions  of  Mach  number,  total  angle  of  attack,  and  re- 
duced azimuth  angle  of  the  angle  of  attack  plane,  0'. 


AERODYNAMIC  ORIENTATION  ANCLES 

In  this  section  the  aerodynamic  orientation  angles  required  for  the 
calculation  of  the  force  and  moment  coefficients  given  in  Equation  40  are 
defined.  The  angle  of  attack,  sideslip,  and  total  angle  of  attack  are  shown 
in  Figure  2.  It  should  be  noted  that,  for  a given  orientation  of  the  body 
with  respect  to  the  velocity  vector,  the  values  of  a and  p depend  upon 
whether  the  body  axes  are  rolling  or  nonrolling.  The  value  of  a,  however, 
is  independent  of  the  choice  of  body  axes.  Also  note  that  and  a always 

lie  in  the  x - z.  plane. 

A A 


a = Arctan 


180°  > a > - 180° 


Range  of  a 


Sign  of  Sign  of 


90°  > a > o° 
180°  > a > 90° 

0°  > a > -90° 
-90°  > a > -180° 


P = Arctan 


180°>p>-180°  (42) 

with  the  same  sign  conventions 
and  quadrants  as  a 


a = Arctan 


,/v  * V * ‘ 

I ^BA  + ^BA 


V 


X 


BA 


180°  > a >0° 


(43) 


Range  of  a 


90°  > a > 0° 
180°  > a > 90° 
a = 0 


Sign  of 
V 

*BA 

Sign  of 
BA 

Sign  of 
V 

^BA 

+ 

^0 

+ 0 

- 

f 0 

t 0 

0 

0 

0 

The  roll  angle  0,  as  is  the  case  with  a and  p depends  upon  whether 
or  not  the  body  frame  is  rolling 


0 = Arctan 


360°  > 0 >0° 


(44) 


Range  of  0 


Sign  of  Sign  of 

V V 

^BA  ^BA 


90°  > 0 > 0°  + 

180°  > 0 > 90°  + 

270°  > 0 > 180° 

360°  > 0 > 270° 

0=0  0 


+ 

0 
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Initially,  the  rolling  and  nonrolling  body  frames  are  coincident. 

If  it  is  desired  to  compute  in  a nonrolling  frame  when  there  is  a con- 
figurational asymmetry  reflected  in  a body-fixed  force  or  moment,  then 
the  precise  roll  position  of  the  actual  body  must  be  computed.  The 
differential  equation  describing  the  roll  relationship  between  the  AB 
(rolling  or  nonrolling  computational  frame)  and  the  body-fixed  frame 
is  Equation  35.  reproduced  here: 

(p'  = (1  - k)p 

However,  this  procedure  is  not  particularly  recommended  since 
the  point  of  introducing  the  nonrolling  computational  frame  is  to  make 
it  unnecessary  to  calculate  the  roll  angle  of  the  body-fixed  frame  ac- 
curately: if  an  asymmetry  must  be  introduced  in  the  body-fixed  frame, 
its  position  must  be  accurately  known.  In  any  case,  whether  the  roll 
position  is  known  accurately  or  not.  there  exists  a relationship  among 
the  roll  angles:  the  azimuthal  position  of  the  angle  of  attack  plane  in 

the  body-fixed  frame  is  equal  to  azimuthal  position  of  the  angle  of 
attack  plane  in  the  rolling  or  nonrolling  computational  frame  plus  the 
roll  angle  between  the  rolling  or  nonrolling  computational  frame  and 
the  rolling,  body-fixed  frame: 

(p  = 0 + cp'  (45) 

Obviously,  when  k = 1,  q)'  = o.  Therefore  cp'  = o and  the  AB  frame 
is  identical  with  the  body-fixed  frame. 

As  an  example  of  the  use  of  roll  symmetry  in  the  calculation  of 
aerodynamic  coefficients . consider  a cruciform  finned  body . For  such 
a body,  the  coefficients  which  are  functions  of  roll  angle  (C^,  c^. 

C , C , and  C ) are  input  from  zero  to  90°  only,  since  these  coefficients 

repeat  their  values  over  each  90°  of  roll.  From  (the  true  azimuthal 
angle  of  the  angle  of  attack  plane  in  the  rolling,  body-fixed  axes)  cp, 
we  calculate  a ficticious  roll  angle  0'  (the  reduced  roll  angle)  . 


0' 

= «P 

90° 

> 

<P 

> 

0° 

0' 

= cp 

-90° 

180° 

> 

<P 

> 

90° 

0' 

= cp 

-180° 

270° 

> 
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> 

180' 

0' 

(p 

-210° 

360° 

> 

9 

> 

210' 

I 


r 


1 


If  the  roll  dependent  coefficients  repeat  every  180°,  the  reduced 
roll  angle  is 


(6'  = (p  180°  > 9 > 0° 

(i'  = cp  -180°  360°  > <p  > 180° 


Provision  is  made  for  a repeat  angle  of  60°,  typical  of  a six- 
finned  projectile,  thus: 


0' 

= 9 

0 

o 

CO 

> 

9 

> 

0° 

= 9 

O 

o 

CO 

1 

120° 

> 

9 

> 

60° 

ci' 

= 9 

-120° 

180° 

> 

9 

> 

o 

CM 

d' 

= <p 

-180° 

240° 

> 

9 

> 

180' 

= 9 

-240° 

300° 

> 

9 

> 

240' 

0' 

= cp 

-300° 

360° 

> 

9 

> 

300' 

All  other  cases  require  that  the  roll  dependence  be  supplied 
from  0°  to  360°,  and  therefore 

0'  = cp 

for  all  values  of  (p . 

In  all  cases  the  table  look-up  and  interpolation  is  performed 
with  d'as  the  independent  variable  for  roll  angle  dependence. 

ROUTINE  INTERMEDIATE  VARIABLE  CALCULATIONS 

The  operation  of  the  program  requires  regular  calculations  of 
dynamic  pressure,  Mach  number,  velocity  of  the  air  frame  with  respect 
to  the  local  air  mass,  and  other  quantities  from  the  fundamental  vari- 
ables X_,  Y Z their  derivatives,  and  the  Matrix  [8„]  . We  define 
W H H B 

a local  axis  system  see  Figure  3,  such  that  its  origin 

is  at  the  burnout  center  of  mass  of  the  body,  its  x axis  pointing 

Lb 

North;  its  z axis  points  East,  and  its  y axis  points  directly  away 
L£  Lb 

from  the  center  of  the  spherical  Earth.  Its  location  on  the  Earth  is 

defined  by  its  north  latitude,  \j/  and  its  east  longitude  t . 

H K 
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The  components  of  the  velocity  of  the  body  with  respect  to  the 
Farth.  in  components  along  the  local  axes  (subscripted  LE)  are 


- — 

1 

~ — , 

V' 

X 

o 

LF 

1 

V 

V 

- 

'^R 

- 

0 j 

^LE 
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\' 
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CO  J y ^ + z * ' 
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E V ^R  R 1 

LE 

^ i 
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_ _l 

(46) 


where 


cos  (v}/j^) 


-sin(\!;j^)  cos(tj^) 


-sin  (>];„)  sin(t^) 
K H 


sin 

K 


cos(ij/„)  cos(t„) 
K K 


cos(vl/j^)  sin(tp) 


-sin  (t^) 

H 


COS(Tp) 


and  CO  is  the  spin  of  the  Earth  and  y and  z are  the  equatorial  com- 
E ^ ^ . 

ponents  of  the  position  of  the  burnout  center  of  mass  in  inertial  coor- 
dinates . 


The  angles  i:  and  are  simply  related  to  the  components  of 

R R 


the  position  vector; 


- ARCTAN  360°  > tp  > 0° 


R 


(47) 


The  usual  quadrant  assignment  for  is  followed  as  a function  of  the 

r ■ 

R' 


signs  of  and  Y^: 


Range  of  t 


R 


Sign  of  Z 


R 


Sign  of  Y 


R 


90°  > > 0° 

R 


180°  > > 90° 

H 


+ 

+ 


270°  > > 180° 

R 


360°  > > 270° 

R 
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and  where  = y„  = o,  t„  is  arbitrarily  defined  to  be  zero. 
R H K 


Similarly 


«j/^  = ARCTAN 
R 


I , ,1  180”  > ♦ > - 180" 

Pr  * “b  / " " 


The  quadrant  assignments  for  \|/j^  are  somewhat  simpler  as  the  positive 
radical  is  always  assumed: 


Range  of  vl/ 

H 


Sign  of  X 


90°  > > 0°  + 

H 

0°  > \!/„  > - 90° 

K 

Since  the  missile  will  not  ever  fly  substantially  below  the  nominal 
surface  of  the  Earth,  no  value  of  need  be  defined  for  = y^  = z^^  = o . 


The  components  of  the  velocity  of  the  vehicle  with  respect  to  the 
local  air  mass  depends,  of  course,  on  the  winds.  The  wind  is  described 
by  a constant  speed  and  azimuth  angle  or  by  a table  of  speeds  and  a table 
of  azimuth  angles,  each  as  a function  of  altitude.  The  wind  is  assumed 
to  be  horizontal  and  the  air  mass  is  assumed  to  move  rigidly  with  the 
Earth,  except  for  this  wind.  Therefore 


V + V cos  (A  ) 
^LE  ^ 


V + V sin  (A  ) 
z,  ^ w w 
LE 


where  A is  the  azimuth  angle  of  the  source  of  the  wind,  measured  posi- 
w 

tively  from  the  North  toward  East. 
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The  components  of  this  air  velocity  must  be  resolved  in  the  rolling 
or  nonrolling  computational  frame.  Symbolically,  this  is 


V 

BA 

""la 

V 

^BA 

V 

^LA 

V 

V 

^BA 

^LA 

— 

_ 

where  subscript  T denotes  the  transpose  of  a'  matrix . These  rotation 
matrices  are  unitary,  so  that  their  transposes  are  their  inverses. 

The  total  velocity  of  the  body  with  respect  to  the  local  air  mass  is 


V = ,1  V * +V  2 ' (51) 

U BA  ^ BA  BA 

and  the  Mach  number  is  given  by 
Va 

M = — (52) 

c 


where  the  speed  of  sound,  c,  is  calculated  from  the  altitude  by  means  of 
the  standard  atmosphere  of  Reference  3 or  an  atmosphere  table  supplied 
by  the  user . 

The  altitude  of  the  vehicle  is 


from  the  assumption  of  a spherical  Earth,  and  the  dynamic  pressure  is 
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The  atmospheric  density  is  calculated  by 
p = (D  + 1) 

where  p^  is  the  atmospheric  density  from  either  the  standard  atmosphere 
or  user  supplied  tables,  and  D is  a density  deviation  table,  versus  altitude, 
supplied  by  the  user. 


INITIAL  CONDITIONS 


The  user  of  a trajectory  program  wishes  to  give  the  initial  condi- 
tions as  measurements  in  an  Earth-fixed  frame  of  reference  at  the  launch 
point,  while  the  program  computes  in  either  inertial  axes  or  body-fixed 

axes.  Therefore  X_,  Y_ , Z their  derivatives,  and  [C  ] must  be  ex- 
R R K ^ D 

pressed  as  functions  of  h,  V^,  A^,  y^,  v|/g,  t^,  0,  a,  d,  p,  q,  andr. 

Since  the  Earth  and  inertial  frames  are  initially  coincident  (Fig  6) 
the  initial  values  of  inertial  latitude  and  longitude  are  equal  to  the  Earth 
frame  values  input.  Inverting  Equations  47,  48,  and  53  allows  us  to  solve 
for  the  initial  coordinates  in  terms  of  the  altitude,  radius  of  the  Earth, 
latitude,  and  longitude. 


The  components  of  , the  velocity  of  the  vehicle  with  respect  to 

the  Earth,  shown  in  Figure  4 resolved  in  the  local  frame,  can  be  calculated 
by  inverting  Equations  56,  57,  and  58  simultaneously: 

360°  > Ag  > 0°  (56) 

where  the  quadrant  assignments  for  A are  again  the  traditional  ones, 

E 


Tv 


= ARCTAN 
E 


LE 


V 


LE 


Range  of 


Sign  of  V 

^LE 


Sign  of  V 

LE 


90°  > A_  > 0° 
E 

+ 

180°  > A_  > 90° 

E 

+ 

270°  > A^  > 180° 
E 

- 

360°  > A„  > 270° 
E 

- 

+ 


+ 
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= 0.0. 


and  the  special  value  of  0.0  is  assigned  to  A when  V = V 

\e  ^LE 


Ye  is 


Ve  = iv  ^ ' 

\ ^LE  ^LE 


- 90°  > > 90° 

— * E — 


where  has  the  sign  of  V and  is  zero  when  V is  zero.  is 

^LE  ^LE 


V^  = \ V2  +V2  +y  2 
^ \ ^LE  ^LE  ^LE 


Once  we  have  obtained  the  components  of  V in  the  local  frame,  the 

£ 

components  in  the  inertial  frame  can  be  calculated  by  inverting  Equation  46, 


The  azimuth  and  elevation  angles  of  the  vehicle's  velocity  vector 

with  respect  to  the  local  air  mass  are  shown  in  Figure  4 and  are  computed 

from  the  local  frame  component  of  the  vehicle's  velocity  with  respect  to 

the  Earth,  V , V , and  V . To  obtain  the  velocity  relative  to  the 
X V z 

LE  ^LE  LE 

local  air  mass,  we  use  Equation  49, solving  for  V , V , and  V 


which  can  then  be  used  for  computing 


\a  ^LA 


A,  = ARCTAN  — 
A V 


360°  > A,,  >0° 
A — 


Ya  “ ARCTAN  1 y 2 + V 2 


90°  >Y^  > + 90° 


^LA  ^LA 


The  sign  and  quadrant  relations  for  A^  and  y^  are  the  same,  respectively, 

as  those  of  A and  Yt:..  defined  for  Equations  56  and  57. 

£ £ 
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The  initial  elements  of  the  rotation  matrix  [6„]  are  also  obtained 

£5 

in  an  inverted  manner.  The  matrix  [^^1  is  the  transpose  of  the  matrix 
product  of  a sequence  of  rotations  transforming  an  axis  system  initially 
coincident  with  the  inertial  axes  to  be  coincident  with  the  rotating  or  non- 
rotating computational  frame,  which  is  initially  coincident  with  the  body- 
fixed  frame.  The  sequence  of  rotations  is 


amount 

about  the 

^R 

axis. 

producing  a frame 

^R*  ^R-  "r- 

(Fig  3) 

amount 

-^R 

about  the 

^R 

axis. 

producing  a frame 

^LE’  ^LE'  ^LE' 

(Fig  3) 

amount 

about  the 

^LE 

axis, 

producing  a frame 

x ' y,'  z ' 

LE  '^LE  LE 

(Fig  4) 

amount 

about  the 

^LE 

axis. 

producing  a frame 

LF’  ^LE’  LE' 

(Fig  4 
6 5a) 

amount 

6-90° 

about  the 

axis. 

producing  a frame 

\e-  ^LE’  LE' 

(Fig  5b) 

amount 

a 

about  the 

V 

^LE 

axis. 

producing  a frame 

^A'  ^A-  "a  ' 

(Fig  5b) 

amount 

d 

about  the 

^A 

axis. 

producing  a frame 

^BA’  ^BA’  ^BA’ 

(Fig  2) 

A positive  rotation  about  the  X-axis  of  amount  6 is  represented  by 
the  rotation  matrix 


[i 


loo 
o cos (6)  sin (6) 

o -sin (6)  cos (8) 


and.  similarly,  positive  rotations  about  the  Y-axis  and  7.-axis  are  repre- 
sented by 


[e  ] 
y 


cos (6)  o 

o 1 

sin  (6)  o 


-sin  (6) 
o 

cos (6) 


cos (6)  sin (6)  o 

[t  1 = -sin (6)  cos (8)  o 
z 


Applying  these  representations  to  the  sequence  of  rotations  described 
above,  we  have  for  the  initial  value  of  the  matrix  transforming  vectors 
from  the  inertial  frame  to  the  rotating  or  nonrotating  computational  frame 


1 o o cos  (a)  o -sin  (a) 

T 

[Sg]  = o cos(ci)  sin(eS)  o 1 o 

o -sin(«S)  cos(0)  sin  (a)  o cos  (a) 


0 

0 

cos(y^) 

sin  (Y^) 

0 

sin  (B) 

-cos (0) 

-sin  (y  . ) 
A 

cos  (Y^) 

0 

cos (0) 

sin  (0) 

0 

0 

1 

0 

sin  (A^ ) 
A 

! [cosCxj/  ) 
1 ^ 

-sin(i!-'j^) 

0 

jsin(vl/g)  cos(v|/g) 


-sin(A  ) o cos  (A  ) 

A A 


o cos(tg)  sin(tg) 


-sin(tg)  cos(Tg) 


Given  the  matrix  transformations  [C„]  and  [5_  ] , the  components 

ID  L 

of  velocity  of  the  vehicle  relative  to  the  local  air  are  calculated  by  Equation  50. 


.•1^  . 


This  then  allows  the  computation  of  the  angles  of  attack  and  sideslip  il- 
lustrated in  Figure  2. 
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a = ARCTAN  (V  /V  ) 

^BA  ^BA 

and 

P = ARCTAN  (V  /V  ) 

^BA  *BA 

SPECIAL  CALCULATIONS  FOR  OUTPUT 

These  calculations  are  made  to  either  convert  the  state  variables 
used  for  the  integration  procedure  into  variables  convenient  for  a user  at 
each  printed  interval  or  for  special  versions  of  the  program.  First  we 
will  consider  the  former . 

The  longitudinal  and  latitudinal  ranges  are  the  great  circle  arc 
length  on  the  surface  of  the  Earth  from  zero  degrees  longitude  and  latitude 
to  the  present  position  of  the  intersection  of  the  radius  vector  of  the  pro- 
jectile with  the  Earth's  surface. 


' "e'e 

(62) 

(63) 

where 


= ARCTAN  (Z^/Y^)  360°  > t > 0° 

b EE  — E 


(64) 


and 


= ARCTAN 


180°  > v|/^  > -180° 
— E 

(65) 


The  sign  conventions  for  the  Earth's  longitude  and  latitude  are  the 
same  as  those  for  Equations  47  and  48.  The  longitudinal  range,  is  always 
positive  and  the  latitudinal  range,  R , has  the  sign  of  vj/^. 
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The  Earth-fixed  coordinates  of  position  of  the  projectile  are  obtained 
from  the  inertial  components  of  the  position  vector  by 


s ~ 

^E 

MV 

"r 

^E 

_ _J 



(66) 


The  Earth  and  inertial  frames  are  initially  coincident  and  separate 
from  each  other  by  the  spin  of  the  Earth,  which  is  about  the  Earth's  X-axis 
(Fig  6) . Therefore,  the  rotation  matrix 


[e. 


cos(co„t)  sin(co_t) 


-sin  (co  t)  cos  (03  t) 
E E 


(67) 


is  defined  for  any  time  of  flight. 

The  time  rate  of  change  of  the  aerodynamic  roll  angle,  «i,  can  be 
useful  when  the  computational  frame  rolls  with  the  body . A small  oscil- 
lating value  of  h,  averaging  zero,  indicates  that  the  projectile  is  in  "lunar 
motion"  or  is  "locked  in,"  a condition  which  indicates  yaw  will  grow. 
Therefore  we  compute  d for  printing . 


Tan(0)  = 


from  Equation  44 


i 


f 


Taking  a derivative  with  respect  to  time. 


V 


- V 


d = cos*  («5) 


BA  •'BA 


V 

^BA  ^BA 


V * 
^BA 
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but 
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therefore 


1 — - r V 
m X, 


1—  + q V ] 


Several  of  the  alternate  versions  of  the  program  require  special 
calculations.  PANOL  3 and  PANOL  5 both  allow  the  center  of  mass  to  lie 
off  the  centerline  of  the  computational  frame  and  allow  the  moment  of  in- 
ertia tensor  to  be  non-diagonal.  Both  versions  of  the  program  are  valid 
only  after  burnout,  when  the  thrust  and  jet  damping  integrals  are  zero 
as  is  the  gravity  moment.  (See  Equation  7.)  The  gravity  moment  induced 
by  the  off-axis  position  of  the  center  of  mass  is  ignored  as  it  averages  to 
zero  due  to  spin.  For  these  conditions  we  can  rewrite  the  rotational  equa- 
tion of  motion  in  tensor-vector  form , as 


where  H is  the  angular  momentum  vector,  the  aerodynamic  torques, 
and  the  aerodynamic  forces;  is  in  this  case  the  position  of  the  center 

of  mass  with  respect  to  the  nominal  burnout  center  of  mass  position  (the 
computational  frame  origin)  . 

If  we  fix  the  coordinate  system  to  the  projectile,  the  time  rate-of- 
change  operator  in  inertial  space  must  be  transformed: 


+ X 


but  the  angular  momentum  becomes  simply 


H = [I]  Q 


where  [I]  is  the  moment  of  inertia  tensor  and  is  constant  in  body-fixed 
coordinates . 


Therefore 


and 


where 


^ -V 

+ F,  xr  +rJxH 
A c 


d_ 

dt 


n = [I] 
BA 


^ ^ 
{r.+F^xr  +P.xH} 
A A c 


(69) 


GUIDANCE  AND  CONTROL  CALCUU^TIONS 

The  first  requirement  for  a controlled  vehicle  is  the  definition  of  a 
target.  Added  to  this  program  is  a rectangular  target  in  a vertical  position. 
The  target  is  of  variable  height  and  length.  It  is  divided  into  a number  of 
subrectangles,  each  of  which  has  an  intrinsic  brightness  defined  by  the 
user.  The  range,  deflection,  and  orientation  of  the  target  can  be  input 
as  tables  versus  time  and  are  linearly  interpolated  for  intermediate  values 
of  time . 

Examining  Figure  8,  we  see  the  target  at  longitude  and  latitude  on 

the  surface  of  the  Earth.  It  has  a length  (EL)  and  a height  (HITE)  and 

its  unit  normal  (n)  is  rotated  by  angle  AZM  from  north  toward  west.  In 
target  coordinates,  the  center  of  a subrectangle  (dA)  is  at  Yj^qq  . 

z . In  Earth  coordinates,  a vector  from  the  origin  to  a given  element 
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r 


of  the  target  (dA)  is 


(Rp,  + Vloc'  ^LOC  ^ii^^AZMIcosCvi/^) 

^dA  ^ ^LCC^  cos(v|/^)cos(t^)  - {cos  (AZM)  sin  (x^J  + 

sin  (AZM)cos  (t„)sin  (»l/_) } i 

T T I 

(Rg  + VLOc^  cos(vt(T)sin  (17)+  Zgoc  {cosfAZMicosft^)  - i 
sin  (AZM)  sin  (t^)  sin  (vl/^) } 

The  position  of  the  sensor  in  the  missile  coordinates  is  at  ^g-Yg-^g 

and  since  for  controlled  vehicles,  the  Earth  is  not  allowed  to  rotate,  the 

matrix  [^^1  transforms  vectors  from  body  frame  effectively  to  Earth  frame. 

Since  R^  is  the  radius  of  the  Earth  and  h the  altitude  of  the  vehicle,  the 
E 

sensor,  in  Earth-coordinates  is  at 


1 

— 

(Rp,  + h)  sin(\|/j^) 

X 

S 

R = 
s 

(Rp,  + h)  cos(\}/p^)cos(Xpj) 

^s 

(Rp,  + h)  cos  (\|(p^)sin  (Xp^) 

2 

S 

Therefore,  the  vector  from  a target  element  to  the  sensor  is 


R = R - R^. 
s dA 


A unit  normal  from  any  area  element  of  the  target  is 


cos  (AZM)  cos 

n=  sin  (AZM)  sin  (X7)  - cos  (AZM)  cos  (x^)  sin  (if/^)  (72) 


sin  (AZM)  cos  (x^)  - cos  (aZM)  sin  (x^)  sin  (vjj^) 


The  visible  or  projected  area  (of  the  differential  area  dA)  in  the 
direction  of  R is 

. R • n 

which  we  multiply  by  the  intrinsic  brightness  of  the  differential  area  and 
divide  by  the  square  of  the  distance  between  the  differential  area  and 
the  sensor . 

The  sensor  is  misaligned  with  the  projectile  X-axis  by  angle  YAW 

in  a plane  rotated  by  angle  ROL  from  the  - y^  plane.  A matrix  [C^] 

produced  by  the  rotation  of  amount  ROL  about  the  X -axis  followed  by  a 

rotation  of  amount  YAW  about  eh  Z ' axis  describes  this  misalignment. 

B 


cos  (YAW) 


sin  (YAW)  o 


-sin  (YAW)  cos  (YAW) 


cos  (ROL)  sin  (ROL) 

-sin  (ROL)  cos  (ROL) 


cos  (YAW)  sin  (YAW)  cos  (ROL)  sin  (YAW)  sin  (ROL) 

-sin  (YAW)  cos  (YAW)  cos  (ROL)  cos  (YAW)  sin  (ROL) 

o -sin  (ROL)  cos  (ROL) 

(74) 


Therefore,  L (a  vector  along  the  line  of  sight)  is,  in  body  frame, 

I 1 n , cos  (YAW) 


L[£J 


1 

O !-• 

. 1 

1 

r 

o 

1 

sin  (YAW)cos  (ROL) 
sin  (YAW)  sin  (ROL) 
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and  in  Earth  frame. 


T 


▼ 


cos  (YAW) 

sin  (YAW)  cos  (ROD 

sin  (YAW)  sin  (ROD 


1 

I 

I 

i 
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Therefore,  the  angle  between  the  line  of  sight  of  the  seeker  and  the 
line  from  the  seeker  to  the  target  area  element  is  ALEPH. 


ALEPH  E ARCCOS 


- L ^ 

I R| 


(75) 


ALEPH  may  or  may  not  be  within  the  acceptance  angle  of  the  sensor  (ACANG)  . 
If  not,  the  element  of  target  is  ignored:  if  within  the  acceptance  angle,  the 
above-mentioned  projected  area  (A  ) is  multiplied  by  the  intrinsic  bright- 
ness and  divided  by  the  square  of  the  distance  to  the  area  element  to  form 
a weighting  factor  (W..)  and  the  vector  R is  resolved  in  sensor  oriented 
components.  If  we  examine  Figure  9 we  see  a single  plane  (x-sensor, 
y-sensor)  of  the  sensor-target  relationship. 


The  y-sensor  component  of  the  image  is  a position  y,  which  is  pro- 
portional to  the  y^  component  of  the  vector  R 


y = KR, 


We  wish  to  scale  the  image  of  the  target  such  that  when  the  i,j^^ 

target  element  is  at  the  limit  of  the  acceptance  angle  its  image  element  is 

at  y..  = 1 . 
ij 


TAN  (ALEPH)  = 


R 


^s/r^^  = TAN  (ACANG) 


since  .. 

y^^  = 1 = K R^  TAN  (ACANG) 
s 

therefore 


TAN  (ACANG) 

X 

s 
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and 


TAN  (ACANG 

X 

s 


The  centroid  is  then  approximated  by 


and 


SUM 


T 

i] 


W.. 

13 


YCNTR  = I y 

i.j 


3V../ 

ij/ 


SUM 


A similar  expression  can  be  derived  for  the  z - direction.  The 
total  strength  of  the  signal  to  the  sensor  is  SUM  and  the  polar  coordinates 
of  the  centroid  of  the  image  brightness  can  be  calculated  from  the  Cartesian 
coordinates  in  the  usual  way. 

An  empirical  formula  for  the  azimuthal  and  radial  error  in  target 
centroid  position  due  to  vehicle  roll  rate  interacting  with  a moving  reticle 
is  included . 

The  differential  equations  of  motion  of  the  gyroscope  were  developed 
by  Dr . Gene  Hopkins  of  Texas  Instruments  for  the  Tube  Launched  Guided 
Projectile.  Much  of  the  description  of  the  gyro  model  is  from  Reference  4. 

The  gyro  comprises  a spherical  gas-bearing  torquers  and  a gimbal- 
angle  readout  system.  Gyro  motion  is  described  with  respect  to  a vehicle 
body-fized  coordinate  system  by  the  angles  0 (produced  by  a positive 

rotation  about  the  vehicle  y-axis  and  \j/^  (produced  by  a positive  rotation 
about  the  [intermediate]  z-axis)  . In  matrix  notation,  the  rotations  take 
the  form 

0 -sin(0  ) 
g 

1 0 

0 cos (0  ) 

g 


cos(0  ) 
g 


sin  (0 


g 
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and 


[T, 


cos  (v[/ 


g 


sin  (vl/^) 


-sin('i'  ) cos(v(/  ) o 
g g 


The  total  ti  ansformation  is 

cos  (i}(  )cos(6 


[T  ] [T  ] 
vj/  6 


g'  g 
sin(v|/  )cos(0 


sin  (0  ) 
g 


sin  (\j>  ) -cos(\)i  ) sin(0  ) 


g 


g 


cos(\|/  ) sin(0  ) sin(v{/  ) 
g g g 

o cos (0  ) 

g 


(76) 


This  vehicle  spin  vector  as  expressed  in  body  frame  is 

(o  = pi  + qj  + rk 

and  the  total  spin  of  the  intermediate  frame  (x',  y' . z')  is 

(Oi  = [Tg]  m + 0g  j 

The  total  spin  of  the  gyro  frame  (that  is,  gyro  spin  less  the  axial 
component  v/ith  respect  to  the  body)  with  respect  to  inertial  space  is 

m = [T  ] coi  + vl/  k 
g '1'  g 

In  gyro  frame 

. ^ A A 

CO  =p  i =qj  +r  k 
g g g g g g g 

Solving  for  components,  we  get 

p =pcos(vl/  )cos(0  ) + (q  + 0 )sin(vl/  ) - r cos  (v]/  )sin(0  ) 
g g g g g S g 

q = - p sin  (vie  )cos(0  ) + (q  + 0 )cos(vl/  ) + r sin  (\{/  )sin(0  ) 
g g g g g g g 


r = p sin(0  ) + r cos(0  ) + \]/ 
g g g g 
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Substituting  the  second  equation  into  the  first  and  eliminating  0^ 

we  get 

p = [p  cos(9  ) - r sin  (6  ) + q sin  (\)/  ) ] sec  (\|/  ) 

^g  g g g g g 

From  the  second  equation  we  can  then  find 

e = [q^  + p sin  cos  (0^)  - r sin  (»|>  ) sin  (0  ) ] sec  - q 


g g g 

and  from  the  third  equation 


vj;  = r - p sin  (0  ) - r cos  (0  ) 
g g g g 


These  show  the  exact  kinematic  description  of  the  gyro  gimbal 
motion.  However,  the  angular  momentum  for  a symmetric  rotor  is 

H=I  (p+co)i+I  qj+I 
g s'  g y^  ^g  'g  y^.  g g 

where  I is  the  axial  moment  of  inertia  and  I the  transverse  moment 
X y_ 

r r . . 

of  inertia  of  the  rotor.  We  can  write  the  equation  of  motion  under  the  in- 
fluence of  applied  torques,  M,as 

dit 


and  then  by  components 
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Since  gyro  rotor  speed  is  considered  constant,  an  expression  for 
CO  is  not  needed.  The  moment  of  inertia  ratio  of  the  rotor  is  chosen  such 

that  the  nutation  frequency  of  the  gyro  is  much  above  the  spin  frequency 
of  the  gyro.  Therefore,  the  only  mode  of  the  gyro  we  need  to  analyse  is 
the  precession,  for  the.qther  components  will  be  outside  the  band  width 
of  the  projectile  and  control  system.  These  higher  frequencies  can  be 
filtered  out  of  the  gyro  pick-off  signal.  We  need  only^^a  low  frequency 
component  of  the  resultant  motion  so  we  assume  q = r = o in  the  moment 

expressions  and  solve  for  q and  r . 

S & 

q=M/[p  (I  -I  )+I  co] 
y/  g y^  x^  s 


r = - M 

g z 


/(Pg  - V * “s' 


r "r  r 


Torque  is  applied  to  the  rotor  by  coils  on  the  outer  frame  of  the 
gyro  so  the  torque  vector  must  be  resolved  into  rotor  coordinates; 


App  = Ty  j + ^2  k 


“ = 1^'  l’'e>  ’App 


M = o by  the  assumption  of  constant  axial  spin 

X 

M = I cos  (\|/  ) + 1 sin  (\!/  ) sin  (0  ) 

y y g z g g 


M = t cos (0  ) 
z z g 


I- 


I 


In  suniH'ary  , we  have 

u = [p  cos(0  ) - r sin  (0  ) + q sin(tl^  )J  seclij/  ) 
g g g g g g 


g z/  xr  ■'r  r 


r„  = M„  / [p  (I  - I ) + I o)  ] 
g y/  g X y x^  s 


r ^ r r 


0 = {q  + [p  cos(0  ) - r sin  (0  )]  sin  (0  )}  sec  (i)/  ) - □ 
g g g g g g 


il/  = r - p sin(0  ) - r cos(0  ) 
g g g g 


for  the  exact  low  frequency  response. 

To  simplify  the  equations  and  avoid  numerical  solution  of  the  algebra, 
we  examine  several  terms.  In  the  equations  for  p^  and  0^,  since  the 

gimbal  angles  are  no  more  than  a few  degrees  and  the  gyro  is  more  nearly 
inertially  still  than  the  vehicle,  q sin  (v|/  ) is  small  compared  with 

r sin(0  ) (on  the  average)  and  very  smafl  compared  with  p cos(0^)  . 

The  p equation  can  then  be  solved  for  p cos  (0  which  can  be  substituted 
g • g 

into  the  0 equation,  leaving 

s 


p a [p  cos(0  ) - r sin  (0  )]/cos(\l/  ) 


0 a q / cos  (vly  ) + p sin  (\1;  ) 
g g/  g g g 


Now  the  gyro  gimbal  angles  can  be  integrated  directly  from  initial 

conditions  and  body  pitch,  yaw,  and  spin  rates.  We  begin  with  the  p 

• g 

equation,  then  q^  and  r^,  finally  using  the  gyro  rates  for  the  0^  and 
equations  and  integrating  by  numerical  means . The  interconnections 


of  the  gyro  loop  may  be  seen  on  Figure  7. 
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The  output  of  the  gyro  loop  in  two  dimensions  is  handed  to  the 
thruster  logic  which  determines  which  of  each  opposed  pair  of  lateral 
thrusters  will  be  commanded  on  if  the  gas  generator  has  been  turned  on 
and  has  not  yet  been  exhausted.  The  valve  diverting  flow  between  the 
opposed  thrusters  is  assumed  to  have  a delay  between  the  command  to 
switch  direction  and  the  execution  of  the  command.  The  flow  actually 
switches  completely  and  instantaneously. 

A display  in  graphic  form  is  provided  at  each  time  printout  is  gen- 
erated. It  displays  a circle  simulating  the  field  of  view  of  the  sensor  in 
body  coordinates,  an  indication  of  the  on  or  off  status  of  each  thruster, 
a mark  indicating  the  centroid  of  the  part  of  the  target  in  view,  the  co- 
ordinates of  the  centroid,  a statement  as  to  whether  control  is  on  or  off. 
and  a statement  as  to  whether  the  signal  from  the  sensor  is  sufficiently 
strong  to  have  the  control  system  locked  on. 
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Fig  4 Velocity  components  in  local  axes  system 
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Fig  7 Guidance  loop  block  diagram 
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j iu«  Tni  Fortran  variable/ svmbol /deter iption 


Fig  8 Target  and  sensor  geometry 
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1 1-7  JO/  / Input  35  to  omit  wind  tables  input 

Input  37  to  include  wind  tables 


8-14  NATMOS/  / Number  of  points  in  nonstandard 
atmosphere  tables;  max  = 30 

15-21  NTHRST/  / Number  of  points  in  thrust  table: 

max  = 30 

22-28  KY  (12)/  / Input  1 to  omit  printing  of  aerodynamic 

tables,  otherwise  input  0 

29-35  KY  (13)/  / Input  1 to  omit  printing  output  formats 

4 and  5 

2 1-80  TYPE/  / Alphameric  title  for  printing  on  output 

Cards  numbered  3 to  16  are  aerodynamic  coefficient  tables  and  can 
be  two  dimensional  or  three  dimensional.  All  tables  are  required  for 
all  versions  of  the  program.  An  example  of  a two  dimensional  table 
is  the  axial  force: 

3-A  1-2  NCXM/  / Number  of  Mach  numbers  for  table; 

max  = 16 

3-4  NCXA/  / Number  of  angles  of  attack  (degrees) 
for  table:  max  = 10 

3-B  1-72  TCXM/  / Table  of  Mach  numbers  for  table, 

E12  format,  6 per  card,  in  ascending  order 

3-C  1-72  TCXA/  / Table  of  angles  of  attack  for  table, 

E12  format,  6 per  card,  in  ascending  order 

3-D-l  1-72  TCX/C  / Table  of  coefficients,  for  lowest  Mach 

number,  one  entry  for  each  angle  of  attack, 
in  order  of  ascending  angle  of  attack, 

6 per  card,  E12  format  (may  be  more  than 
one  card) 
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Card  No. 


Col 


Fortran  variable/ symbol/description 


r 


3-D-2  1-72  TCX/Cj^/  Table  of  coefficients  for  second  Mach 

to  number  as  card  3-D-l  above 

3-D-NCXM  Repeat  until  NCXM  tables  are  supplied 

An  example  of  a three  dimensional  table  is  side  force; 


5-A 

1-2 

NCYM/  / 

Number  of  Mach  numbers  for  C table; 

y 

max  = 6 

3-4 

NCYA/  / 

Number  of  angles  of  attack  for  C table; 
max  =30  ^ 

5-6 

NCYR/  / 

Number  of  roll  angles  for  C table; 
max  = 30  ^ 

5-B 

1-72 

TCYM/  / 

Table  of  Mach  numbers  for  C table 

y 

5-C 

1-72 

TCYA/  / 

Table  of  angles  of  attack  (degrees)  for 
C table 

y 

5-D 

1-72 

TCYR/  / 

Table  of  roll  angles  (degrees)  for  C^  table 

5-E-l 

1-72 

TCY/C  / 

y 

Table  of  C coefficients  for  lowest  Mach 
number  anS  angle  of  attack,  one  entry 
for  each  roll  angle,  in  order  of  ascending 
roll  angle,  6 per  card,  E12  format  (may  be 
more  than  one  card) 

5-E-2 

to 

5-E-NCYA 

1-72 

TCY/C  / 

y 

Table  of  C for  lowest  Mach  number  and 
second  lowest  angle  of  attack  as  card 
5-E-l  above 

Repeat  until  NCYA  tables  are  supplied; 

then  repeat  5-E-l  to  5-E-NCYA  for 
ascending  values  of  Mach  number  until 
NCYM  groups  of  tables  are  supplied 
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The  following  table  gives  the  card  numbers,  coefficient  names, 
dimension  (dim)  and  maximum  number  of  Mach  numbers,  angles  of 
attack  and  roll  angles  allowed  for  each  required  table,  in  order. 


Card 

No. 


Coefficient  name  (symbol) 


Max  no.  Roll 
Dim  Mach  Alpha  angle 


3 

Axial  force  (C  ) 

X 

2 

16 

10 

4 

Spin  induced  axial  force  (C  ) 

2 

6 

5 

5 

Induced  side  force  (C  ) 

3 

6 

10 

30 

y 

6 

Normal  force  (C  > 
z 

3 

16 

10 

15 

7 

Magnus  force  (C  ) 

2 

16 

10 

P 

8 

Trim  force  (C  ) 

2 

6 

5 

^dA 

9 

Induced  roll  moment  (C^) 

3 

6 

10 

30 

10 

Roll  moment  due  to  fin  cant  (C  ) 

2 

6 

10 

‘6 

11 

Roll  damping  moment  (C  ) 

2 

16 

10 

P 

12 

Pitching  moment  (C  ) 
m 

3 

16 

10 

15 

13 

Induced  side  moment  (C  ) 

m 

3 

6 

10 

30 

14 

Pitch  or  yaw  damping  moment 

q 

2 

16 

10 

15 

Magnus  moment  (C  ) 

m 

2 

16 

10 

16 


Trim  moment  (C 


m 


6A 


Card  No. 


Col 


Fortran  variable/ symbol /description 


18 

19 

20 
21 
22 

23 

24 

25 

26 


1 27 


1-2 

3-4 

5-6 

1-72 

1-72 

1-72 

1-72 

1-72 

1-72 

1-72 

1-72 

1-72 

1-72 


NHD/  / Number  of  points  in  atmospheric  density 
deviation  table;  max  = 30 

NHW/  / Number  of  points  in  wind  velocity  table: 

max  = 30.  Not  needed  for  JO  = 35  (card  1) 

NHA/  / Number  of  points  in  wind  azimuth  table; 

max  = 30.  Not  needed  for  JO  = 35  (card  1) 

THD/  / Table  of  altitudes  for  atmospheric  density 

deviation  table,  6E12  format:  max  of  5 cards 

TON/  / Table  of  atmospheric  density  deviations, 

6E12  format:  max  of  5 cards.  Density 
= Density  (1  + deviation) 

THW/  / Table  of  altitudes  for  wind  table  (ft)  , 

6E12  format;  max  of  5 cards 

TWD/  / Table  of  wind  velocities  (ft/sec)  , 6E12 
format;  max  of  5 cards 

THA/  / Altitude  table  for  wind  azimuth  (ft) , 6E12 
format:  max  of  5 cards 

TA/  / Wind  azimuth  table  (deg)  , 6E12  format; 
max  of  5 cards 

ALT/  / Table  of  altitudes  for  non-standard  at- 
mosphere (ft) : omit  if  NATMOS  = 0 on 
card  #1;  6E12  format 

TEMP/  / Table  of  temperatures  for  non-standard 
atmosphere  (Rankine);  6E12  format 

DENN/  / Table  of  densities  for  non-standard  at- 
mosphere (lb-sec^ /ft*  or  slugs/ft^); 

6E12  format 

PRESS/  / Table  of  pressures  for  non-standard 
atmosphere  (Ib/ft*);  6E12  format 
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Card  No.  Col  Fortran  variable/symbol/description 

28  1-72  TIME/  / Table  of  times  for  thrust  table  (seconds)  , 

6E12  format,  omit  if  NTHRST  = 0 on 
card  #1;  max  of  5 cards 

29  1-72  THRUST/  / Table  of  thrusts  for  above  times, 

6E12  format:  same  number  of  cards 
as  TIME  table 

30  1-7  INTT/  / Input  1 for  unpowered  flight;  input  2 

for  powered  flight;  17  format 

8-14  NPUNCH/  / Input  1 to  punch  cards  for  FILM 
read  program;  omit  for  no  cards 
desired 

15-21  NOGEES/  / Input  1 for  no  gravity,  omit  for 
standard  gravity:  17  format 

22-28  KAUF/  / Input  0 for  no  tape  output;  input  1 for 
a tape  of  kinematic  variables  to  start  at 
launch;  input  2 to  start  tape  at  apex; 

17  format 

29-35  IGP/  / Input  1 to  calculate  geographic  position 
on  spherical  Earth;  omit  for  flat  Earth 
approximation:  if  IGP  = 0,  Y (217)  = 

Y (220)  = Y (221)  = 0 

36-42  KY  (15)/  / Input  a positive  number  for  tables  of 
target  range  and  crossrange  to  be 
input  on  cards  42  5 43  in  feet;  input 
zero  for  a single  value  of  each  to  be 
input  on  card  51  in  feet;  input  a 
negative  number  if  latitude  and  longi- 
tude are  in  radians 

30  43-49  NUMPTS/  / Number  of  points  in  target  position 

tables  if  KY  (15)  ^ 0;  maximum  of  30 

50-56  NHT/  / Number  of  vertical  subdivisions  of  the 
target  for  the  brightness  pattern 
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Card  No.  Col  Fortran  variable/symbol/description 

57-63  NLNTH/  / Number  of  horizontal  subdivisions 
of  the  target  for  brightness  pattern 

64-70  KY  (16)/  / Plot  frequency  for  seeker  field  of 

view  and  thruster  logic  state  printed 
on  logical  unit  4;  plots  and  every 
KY  (16) th  printed  step 


31 


10 17  Format 

1-7  KY  (2)/  / Run  number 

8-14  KY  (3)/  / Mode  of  integration:  Input  1 for  Runge- 
Kutta,  input  2 for  Adams-Moulton;  in- 
put 3 for  Runge-Kutta  with  truncation 
error  estimate;  Adams-Moulton  is 
recommended 


15-21  KY  (4)/  / Number  of  differential  equations  for 
integration  in  principal  loop;  always 
input  19 

22-28  KY  (5)/  / Always  omit 


29-35  KY  (6)/  / Always  input  1 


36-42  KY  (7)/  / Input  1 for  including  winds,  omit  for 
no  winds 


43-49  KY  (8)/  / Input  1 to  periodically  correct  the 

numerically  introduced  non-orthogonality 

of  the  [C_]  matrix:  omit  for  no  correction 
13 

50-56  KY  (9)/  / Print  frequency  for  output:  prints  each 
KY  (9)  th  integration  step 

57-63  KY  (10)/  / Print  frequency  of  [C  ] matrix  to 

13 

examine  orthogonality 
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Card  No. 


Col 


Fortran  variable/ symbol /descript  ion 


64-70  KY  (11)/  / Aerodynamic  symmetry  options; 

Input  1 for  90°  symmetry;  Input  2 
for  180°  symmetry;  Input  3 for  360° 
symmetry;  input  4 for  60°  symmetry; 
use  2 for  a body  of  revolution 

All  cards  are  5E  14.6  FORMAT  until  further  notice 

32  1-14  Y (200) /h^/  Terminal  altitude,  where  calculations 

are  to  cease  (feet) 

15-28  Y (201) /t^/  Terminal  time,  when  calculations  are 
to  cease  (seconds) 

29-42  Y (202)/ a^/  Terminal  total  angle  of  attack,  above 

which  calculations  cease  (degrees) 

43-56  Y (203) /Ajj/  Altitude  above  which  atmosphere  is 

assumed  to  have  no  effect  (feet) 

57-70  Y (204) /co  / Rotation  speed  of  the  Earth  (radians/ 

second);  .00007292  is  rotating, 
zero  if  not 

33  1-14  Y (205)/Error/  Input  O.for  integration  interval 

to  be  adjusted  to  limit  relative 
truncation  error;  input  -1.0  to 
limit  absolute  truncation  error 

15-28  Y (206) /ENE/  Input  0.0  to  integrate  with  a fixed 

time  step  equal  to  Y (207);  Input 
a positive  value  to  adjust  the  inte- 
gration step  so  that  the  interval  lies 

, ^ ,„-ENE-3  ,„-ENE 

between  10  and  10  seconds 

29-42  Y (207)/At^/  Initial  interval  of  integration  (seconds) 

43-56  Y (208) /t  / Initial  time  (seconds) 
o 


69 


I 


Card  No. 


Col 


Fortran  variable/ symbol /descript ion 


57-70 

Y 

(209)/K/  Input  0.0  for  computational  frame  to 
be  non-spinning:  Input  1.0  for  com- 
putational frame  to  spin  with  body 

34 

1-14 

Y 

(215) /h/  Initial  altitude  (feet) 

15-28 

Y 

(216) /V  / Initial  speed  with  respect  to  the 
Earth  (ft/ sec) 

29-42 

Y 

(217) /Ag/  Azimuth  of  launch  velocity,  measured 
from  North  toward  East  (degrees) 

43-56 

Y 

(218) /y^/  Elevation  angle  of  launch  velocity 
(degrees) 

57-70 

Y 

(219) /0/  Determines  the  direction  of  the  plane 
formed  by  the  missile  X-axis  and  the 
total  wind  vector  if  the  initial  angle 
of  attack  is  not  zero,  see  Figure  5 
(degrees) 

35 

1-14 

Y 

(220) /t/  Initial  longitude,  if  geographic  position 
output  desired  (degrees) ; otherwise  omit 

15-28 

Y 

(221) /\|//  Initial  latitude,  if  geographic  position 
output  desired  (degrees):  otherwise  omit 

29-42 

Y 

(222) /a/  Initial  total  angle  of  attack  (degrees) 

43-56 

Y 

(223) /d/  Initial  roll  angle  of  plane  of  angle  of 

attack  in  computational  frame  (degrees) 

57-70 

Y 

(224) /p/  Initial  spin  rate  (radians/second) 

36 

1-14 

Y 

(225) /q/  Initial  pitch  rate  (radians/second) 

15-28 

Y 

(226) /r/  Initial  yaw  rate  (radians/second) 

29-42 

Y 

(227) /V  / Wind  speed  (feet/ second) 
w 

I 

( 

i 

[ 

i 

I 

I 

I 


Card  No. 


Fortran  variable/ symbol /description 


43-56  Y (228) /A^/  Wind  azimuth  (degrees) 

57-70  Y (236) /6/  Fin  cant  angle  for  spin  acceleration 

coefficient  (Same  units  as  aerodynamic 
coefficient) 

1-14  Y (230) /I  / Axial  moment  of  inertia  (slug-feet 
squared) 

15-28  Y (231) /I  / Transverse  moment  of  inertia  (slug- 
^ feet  squared) 

29-42  Y (232) /I  / Transverse  moment  of  inertia  (slug- 
feet  squared) 

43-56  Y (233) /m/  Mass  of  vehicle  at  burnout  (slugs) 

57-70  Y (234) /A/  Vehicle  reference  area  (square  feet) 

1-14  Y (235) /d/  Vehicle  reference  length  (or  diameter) 

(feet) 

15-28  Y (237) /6^/  Pitch  aerodynamic  asymmetry  angle 

(Same  units  as  C ) 

n’c  A 
6A 

29-42  Y (238) /6g/  Yaw  aerodynamic  asymmetry  angle 

(Same  units  as  C ) 

A 

6A 

43-56  Y (257) /x  / Axial  shift  of  burnout  center  of  mass 

due  to  vehicle  unbalance  (inches) 

57-70  Y (258) /y  / Laterial  shift  of  burnout  center  of 

mass  due  to  vehicle  unbalance  (inches) 


Y (259) /X/  Azimuthal  orientation  of  lateral  center 

of  mass  position  due  to  vehicle  unbalance 
(degrees) 


. 


Card  No. 


Col 


Fortran  variable/ symbol /description 


15-28  Y (261) /Hj/  First  Euler  angle  in  transformation 

from  computational  frame  to  principal 

axis  frame;  rotation  about  the  Z , „ 

^ AB 

axis  (degrees) 

29-42  Y (262) /H2/  Second  Euler  angle:  about  the 

Y'  axis  (degrees) 

43-56  Y (263) /H3/  Third  Euler  angle;  about  the  X” 

axis  (degrees) 

Cards  40-44  are  input  for  powered  flight  only: 
INTT  on  card  30-2 

40  1-14  Y (210) /ti/  Time  when  ignition  occurs  (seconds) 

15-28  Y (211) /t2/  Time  when  burnout  occurs  (seconds) 

29-42  Y (212) /m^/  Initial  mass  of  vehicle  (slugs) 

43-56  Y (213) /I  / Initial  axial  moment  of  inertia  of 

0 vehicle  (slug-ft^ ) 

57-70  Y (214) /I  / Initial  transverse  moment  of  inertia 
^o  of  vehicle  (slug-ft* ) 

41  1-14  Y (240)/ 61/  Thrust  malalignment  in  yaw  plane 

(radians) 

15-28  Y (241) /g2/  Thrust  malalignment  in  pitch  plane 

(radians) 

29-42  Y (242) /X  / Axial  location  of  nozzle  exit  in  body 

for  thrust  malalignment  (feet) 

43-56  Y (243) /Y  / y-lateral  location  of  nozzle  for  thrust 
0 

malalignment  (feet) 

57-70  Y (244) /Z  / z-lateral  location  of  nozzle  for  thrust 
0 

malalignment  (feet) 
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Card  No. 
42 


Col  Fortran  variable/symbol/description 

1-14  Y (245) /m/  Mass  loss  rate  due  to  propellant 

burning,  normally  negative  (slugs/ sec) 

15-28  Y (246)/m  / Initial  mass  of  propellant  (slugs) 

^o 

29-42  Y (247) /p  / Pressure  of  propellant  gas  at  exit 
of  nozzle  (pounds/ft’) 

43-56  Y (248) /A^/  Nozzle  exit  area  (feet  squared) 

57-70  Y (249) /V  / Velocity  of  propellant  gas  at  exit 

(feet/ second) 


43 


44 


1-14  Y (250) /Xj  / Initial  longitudinal  position  of  front 
° face  of  propellant  grain  with  respect 
to  origin  of  body  frame  (feet) 

15-28  Y (251)/X2  / Initial  longitudinal  position  of  center 
o o!  mass  of  propellant  grain  (feet) 

29-42  Y (252) /Xj  / Initial  longitudinal  position  of  rear 
° face  of  propellant  grain  (feet) 

43-56  Y (253) /X4/  Longitudinal  position  of  nozzle  exit 
plane  for  jet  damping  (feet) 

57-70  Y (254) /ri/  Outer  radius  of  propellant  grain  (ft) 

1-14  Y (255) /rj  / Initial  inner  radius  of  propellant 
° grain  (feet) 

15-28  Y (256) /p  / Density  of  propellant  grain  (slugs/ 

P feet’) 


29-42  Y (239) /X  / Initial  axial  location  of  entire  vehicle 

Q 

o center  of  mass  (feet) 


43-56  Y (260)/  / Specific  impulse  of  propellant  if  mass 

loss  rate  is  to  be  calculated  (Ib-sec/lb) 
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Cards  45  to  52  for  PANOL  5 only 


Use  6E12.0  Format  Until  Further  Notice 
Card  No.  Col  Fortran  variable/ symbol/description 

45  1-72  TARTIM/  / Table  of  times  for  target  location,  six 

per  card,  NUMPTS  points  total,  see 
card  30;  max  of  5 cards  (seconds) 

46  1-72  TARAZM/  / Table  of  azimuth  angles  for  the  normal 

to  the  target  surface  out  of  the  back 
side  (degrees) 

47  1-72  TARPSI/\l/,j,/  Table  of  ranges  to  target  versus  time, 

feet,  if  KY  (15)  on  card  30  is  positive: 
table  of  latitude  angles  versus  time, 
radians,  if  KY  (15)  is  negative 

t 48  1-72  TARTAU/x  / Table  of  crossrange  distance  to  target 

* (feet)  , if  KY  (15)  is  positive;  table 

of  longitude  angles  versus  time  (radians) 
if  KY  (15)  is  negative 

Table  of  brightness  of  the  target  sub- 
rectangles (energy  output  per  square 
foot)  from  left  to  right,  starting  at  the 
of  the  target  and  working  upward  in 
rows.  It  is  suggested  that  a row  at 
top  and  at  bottom  and  a column  at  left 
and  at  right  be  held  at  zero  brightness. 
The  silhouette  of  the  target  can  be 
traced  inside  this  region  and  an  estim- 
ate of  the  location  of  a hit  can  be  made 
by  the  program 

^'eight  of  the  target  (feet) 

' Length  of  the  target  (feet) 

25-36  YAWSNSR/  / Sensor  malalignment  with  the  vehicle 

X-axis  (radians) 
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49 


BRITAR/  / 


50 


1-12  HTTAR/  / 
13-24  ALNTAR/  / 


- . 


I 


Card  No.  Col  Fortran  variable/ symbol/description 

37-48  ROLSNSR/  / Roll  angle  of  alignment  plane  of 

sensor  (radians) 

49-60  ACANG/  / Half-apex-angle  of  the  cone  of  the 

field  of  view  of  the  sensor  (degrees) 

50  61-72  SIGMIN/  / Minimum  received  signal  at  the  seeker 

for  lockon-  (energy  per  square  foot): 
must  contain  allowances  for  size  and 
efficiency  of  seeker  as  well  as  the 
chosen  threshold 

1 

51  1-12  AZMTAR/  / Constant  value  for  azimuth  of  target  I 

if  no  table  input  on  card  46  ' 

13-24  TAUTAR/  / Constant  value  for  crossrange  distance  j 

to  target  if  no  table  input  on  card  48  i 

25-36  PSITAR/  / Constant  value  for  range  to  target  if 
no  table  input  on  card  47 

37-48  XSNSR/  / X position  of  sensor  in  body  coordinates 
(feet) 

49-60  YSNSR/  / Y position  of  sensor  in  body  coordinates 
(feet) 

61-72  ZSNSR/  / Z position  of  sensor  in  body  coordinates 
(feet) 

52  1-12  THRLVL/  / Thrust  of  a single  steering  jet  when 

turned  on  (pounds) 

13-24  TIMON/  / Time  steering  jets  turn  on  (seconds) 

25-36  TDEL/  / Time  delay  in  steering  jet  valves;  time 
from  a command  to  switch  thrust  (eg, 
from  left  to  right)  until  the  command  is 
carried  out  (seconds) 
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r n 


Card  No. 

Col 

Fortran  variable/ symbol/description 

37-48 

TIMOFF/  / Time  steering  jets  shut  off  (seconds) 

49-60 

XTHRUST/  / Longitudinal  position  of  steering  jets 
with  respect  to  the  center  of  mass, 
positive  is  toward  the  nose  (feet) 

Cards  53  to  57  are  for  the  TLGP  version  of 

PANOL 

5 only 

53 

1-12 

HGYRO/  / Time  step  for  gyroscope  integration 
(seconds) 

13-24 

TGYRON/  / Time  gyro  turned  on  (seconds) 

25-36 

TGYROFF/  / Time  gyro  turned  off  (seconds) 

37-48 

GIX/I  / Axial  moment  of  inertia  of  gyro  rotor 
^R  (slug-ft*) 

49-60 

GIY/I  / Transverse  moment  of  inertia  of  gyro 
R rotor  (slug-ft* ) 

61-72 

GlZ/Iy  / Transverse  moment  of  inertia  of  gyro 
R rotor  (slug-ft* ) 

54 

1-12 

GYRMAS/M  / Mass  of  the  gyro  rotor  (slugs) 
R 

13-24 

GYRX/  / X-position  of  gyro  rotor  in  body  frame 
(feet) 

25-36 

GYRY/  / Y-position  of  gyro  rotor  in  body  frame 
(feet) 

37-48 

GYRZ/  / Z'position  of  gyro  rotor  in  body  frame 
(feet) 

49-60 

ALEPH/  / Initial  pitch  angle  of  gyro  rotor  (degrees) 

61-72 

BET/  / Initial  yaw  angle  of  gyro  rotor  (degrees) 

Card  No. 

Col 

Fortran  variable/ symbol /description 

55 

1-12 

GYRDX/  / Misalignment  of  rotor  center  of  mass 
ahead  of  center  of  support  (feet) 

13-24 

GYRDY/  / Misalignment  of  rotor  center  of  mass 

in  Y-direction  from  center  of  support  (ft) 

25-36 

GYRDZ/  / Misalignment  of  rotor  center  of  mass 

in  Z-direction  from  center  of  support  (ft) 

37-48 

GYRWX/  / X-component  of  rotor  initial  spin  rate 
(radians/ second) 

49-60 

GYRWY/  / Y-component  of  rotor  initial  spin  rate 
(radians/second) 

61-72 

GYRWZ/  / Z-component  of  rotor  initial  spin  rate 
(radians/  second) 

56 

1-12 

FKA/K^/  Feedback  loop  preamplifier  gain  in 
auto-pilot 

13-24 

FKT/K^/  Torque  gain  and  effectiveness 

25-36 

FKPG/K  / Gyro  attitude  pickoff  gain 
G 

37-48 

FKS/K^/  Error  signal  amplifier  gain 

49-60 

AVGNGP/  / Mean  of  random  noise  in  gyro  pickoff 

61-72 

SDNGP/  / Standard  deviation  of  random  noise  in 
gyro  pickoff 

57 

1-12 

EPSGL/e  / Maximum  error  signal  (at  which  torques 
become  nonlinear) 

13-24 

EONE/ej/  Maximum  linear  torque  commanded  out 
of  torque  rate  gain 

25-36 

ETWO/e2/  Maximum  torque  commanded  by  torque 
rate  gain 
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APPENDIX  B 

Output  Format,  Titles,  and  Units 
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FORMAT  1 

TIME  H VE  VA 


(seconds)  altitude 
(feet) 


speed  with 
respect  to 
the  Earth 
(feet/ second) 


speed  with 
respect  to 
the  local 
air  mass 
(feet/ second) 


QBAR  GAMMA  E 


dynamic 
pressure 
(pounds/ 
square  ft) 


elevation 
angle  of 
velocity  with 
respect  to  the 
Earth  (degrees) 


ALPHA  BAR  PHI 


FORMAT  2 

TIME 

(seconds) 


total  angle 
of  attack 
(degrees) 


(degrees) 

BETA 

yaw  angle 
(degrees) 


PHI  PRI 
(degrees) 

Q 

pitch  rate 
(radians/ 
second) 
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MA 

Mach  number 
(dimensionless) 


ALPHA 

pitch  angle 
(degrees) 


R 

yaw  rate 

(radians/ 

second) 


J 


FORMAT  3 


f 

i 


f 

i 

I 


TIME 

DELT 

P 

RTAU 

(seconds) 

time  step 
(seconds) 

spin  rate 

(radians/ 

second) 

deflection  to 
the  right  or 
longitude 
(feet  or 
degrees) 

PHIDC  r 

(radians/ 

second) 

FORMAT  4 

TIME 

FTC 

FY 

F^ 

(seconds) 

force  in 

^AB 

direction 

(pounds) 

force  in  the 
Y.„  direction 
(pounds) 

force  in  the 
direction 
(pounds) 

MY 

MZ 

torque  in 

the  Y^  ^ 
AB 

direction 

(pound-ft) 

torque  in 

the  Z . „ 
AB 

direction 

(pound-ft) 

FORMAT  5 


TIME  DDXR  DDYR  DDZR 


(seconds)  accelera-  acceleration 
tion  in  in  y-direction 
earth  fixed  (ft/ second/ 
X-direction  second) 

(ft/ second/ 
second) 


acceleration 

z-direction 

(ft/second/ 

second) 
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I 


RPSl 

range  or 
latitude 
(feet  or 
degrees) 


MX 

torque  in  the 

X , _ direction 
AB 

(pound-ft) 


DXR 

velocity  in 
X-Earth 
direction 
(ft/second) 


J 


FORMAT  5 (Cont'd) 


DYR  DZR 

velocity  in  velocity  in 
Y-Earth  Z-Earth 

direction  direction 
(ft/ second)  (ft/ second) 


FORMAT  6 


TIME 

TH 

WGT 

XXI 

(seconds) 

thrust 

(pounds) 

mass 

(slugs) 

axial  moment 
of  inertia 
(slug-ft- 
squared) 

ZZI 

XC 

transverse 
moment  of 
inertia 
(slug-ft- 
squared) 

center  of 
mass  (feet) 

FORMAT  7 

TIME 

RANGE 

DEFLECTION 

HEIGHT 

(seconds) 

range  or 
latitude 
(ft  or 
degrees) 

deflection  or 
longitude 
(ft  or 
degrees) 

altitude 

(feet) 

FX/M 

FY/M 

FZ/M 

axial 
acceler- 
tion  (ft/ 
sec/sec) 

normal 
acceleration 
(ft/ sec/ sec) 

side 

acceleration 

(ft/sec/sec) 

YYl 

transverse 
moment  of 
inertia 
(slug-ft- 
squared) 
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FORMAT  8 


TIME 

P 

Q 

R 

PDOT 

(seconds) 

spin  rate 

(radians/ 

sec) 

pitch  rate 
(radians/ 
sec) 

yaw  rate 

(radians/ 

sec) 

spin  acceler- 
ation (radians/ 
sec/sec) 

QDOT 

ROOT 

VXEB 

pitch 

accelera- 

tion 

(radians/ 
sec/ sec) 

yaw  accel- 
eration 
(radians/ 
sec/sec) 

body  X com- 
ponent of 
velocity  with 
respect  to 
the  Earth 
(ft/second) 

FORMAT  9 

(PANOL  5 only) 

TIME 

RNGTOTGT 

SEEKR  EPS  Y 

SEEKR  EPS  Z 

THR  1 

(seconds) 

range  to 

target 

(ft) 

seeker  error 
signal  y com- 
ponent 
(degrees) 

seeker  error 
signal  z com- 
ponent 
(degrees) 

steering  jet 
#1  thrust 

THR2 

THR3 

THR4 

LOCKON 

steering  jet 
#2  thrust 

steering  jet 
#3  thrust 

steering  jet 
#4  thrust 

logical  variable 
showing  signal 
strength  limit 
reached 

SIG  STR 

ALEPH 

BET 

signal 
strength 
received 
by  sensor 
(energy/ 
second/ft* ) 

gyro  gymbol 
pitch  angle 
(degrees) 

gyro  gymbol 
yaw  angle 
(degrees) 
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A.  Analytic  Symbols 


LIST  OF  SYMBOLS 


Absolute  acceleration  of  a mass  particle  within  the 
geometric  boundary  of  the  vehicle 

Absolute  acceleration  of  the  origin  of  the  body 
axes  frame 


a ,a  .a 
o o o 
X y z 


Components  of  a^  in  the  body  frame 


Acceleration  of  a mass  particle  within  the  geometric 
boundary  of  the  vehicle  measured  with  respect  to 
the  rolling  body  frame 

Vehicle  reference  area 


Azimuth  of  V.  in  local  frame 
A 


Azimuth  of  V in  local  frame 
£ 


Rocket  nozzle  exit  area 


Azimuth  of  local  wind 


AZMTAR 


Azimuth  of  n,  the  unit  normal  from  the  back  of  the 
target 


Speed  of  sound 


S’^Y’^Z 


Aerodynamic  force  coefficients  defined  in  the 
body  frame 


Aerodynamic  moment  coefficients  defined  in  the 
body  frame 


C,.C  -C 
1 m n 


Aerodynamic  moment  coefficients  defined  in  the 
aerodynamic  data  frame 


x/0(pd/2V^) 


= ac 

" ®^y/a(pd/2V^) 

= ^*"z/a6A 

^ ^*"£/a6 

" ®^e/a(pd/2v^) 

" ^‘^n/a(pd/2V^) 

' ^'=m/3(q<I/2V^)  - 6C 

A 

= 9C 

m/aSA 

Vehicle  reference  diameter 
Atmospheric  density  deviation  function 
Resultant  external  force  acting  on  the  vehicle 
Aerodynamic  forces  defined  in  the  body  frame 
Acceleration  due  to  gravity  at  earth's  surface 

Altitude  of  vehicle  above  earth's  surface 
Moment  of  inertia  tensor  (not  necessarily  diagonal) 
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Axial  moment  of  inertia  of  vehicle  at  burnout 


^ 

i.  j.  k 


m 

m 

F 

m 

F 

m 


Transverse  moment  of  inertia  of  vehicle  about  the 

Y axis  at  burnout 
B 

Transverse  moment  of  inertia  of  vehicle  about  the 

axis  at  burnout 
B 

Unit  vectors  in  body  axes  frame 

Unit  vector  in  the  direction  of  the  line-of-sight 
of  the  sensor 

Body  to  inertial  frame  transfer  matrix 

Inertial  to  earth  frame  transfer  matrix 

Inertial  to  local  frame  transfer  matrix 

Body  to  sensor  frame  transfer  matrix 

Mass  of  vehicle  at  burnout 
Instantaneous  propellant  mass 

Instantaneous  mass  of  propellant  gases  within 
geometric  boundary  of  vehicle 

m + m 

P 


m Propellant  mass  flow  rate,  usually  negative 

M Mach  number 


V , V , V 

^LA  ^LA 

Components  of  (velocity  with  respect  to  the 

local  air  mass)  in  local  frame 

V , V , V 

^LE  ^LE  ^LE 

Components  of  V (velocity  with  respect  to  the 
£. 

Earth)  in  local  frame 

X 

c 

Longitudinal  location  of  instantaneous  center  of 
gravity  of  entire  vehicle  with  respect  to  body  frame 

X , y , z 
e e e 

Location  of  a point  on  the  nozzle  exit  plane  with 
respect  to  the  rolling,  body  frame 

X , y . z 
e e e 

Location  of  centroid  of  nozzle  exit  area  with  respect 
to  the  rolling,  body  frame 

X . y . z 
o o o 

Position  of  center  of  gravity  of  vehicle  measured  in 
the  rolling  body  frame 

Xi  , Xj 

Longitudinal  position  of  front  and  rear  faces  of  pro- 
pellant grain  with  respect  to  origin  of  body  frame 

Xj 

Longitudinal  position  of  instantaneous  center  of 
gravity  of  propellant  grain  with  respect  to  origin 
of  body  frame 

X4 

Longitudinal  position  of  nozzle  exit  plane  with 
respect  to  the  origin  of  the  body  frame 

^A-  ^A 

Aerodynamic  data  frame  axes 

Rolling  or  nonrolling,  body  frame  axes 

^E-  ^E 

Earth  frame  axes 

\e'  ^LE’  ^le 

Local  frame  axes 

90 


Target  frame  axes 


^LOC'^LOC’^OC 


•^R  • ^R-  ^R 


Inertial  frame  axes,  also  position  of  vehicle  center 
of  gravity  with  respect  to  the  inertial  frame 


a 


Angle  of  attack 


P 

Y/ 


6 

6A.  6B 


£i  . 62 


e 


Total  angle  of  attack 
Angle  of  sideslip 

Elevation  angle  of  in  local  frame 

Elevation  angle  of  V in  local  frame 
h 

Average  fin  cant  angle 

Effective  configurational  asymmetry  angle  in 
pitch  and  yaw  planes  of  rolling,  body  frame 

Thrust  misalignment  angles  in  the  yaw  and  pitch 
planes  of  rolling , body  frame 

Orientation  of  plane  of  total  angle  of  attack  (a)  and 
plane  of  elevation  of  (y^) 

Actual  atmospheric  density 

Standard  atmospheric  density 

Density  of  propellant  gases  at  nozzle  exit  plane 

Density  of  propellant  gases  within  geometric 
boundary  of  vehicle 

Density  of  propellant  grain 
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Earth  longitude 


Inertial  longitude 
Target  longitude 

Azimuthal  orientation  of  rolling,  body-fixed  frame 
with  respect  to  -Z^  plane  (angle  of  attack  plane) 

Azimuthal  orientation  of  rolling  body  frame  with 
respect  to  nonrolling  body  frame 

Azimuthal  orientation  of  rolling  or  nonrolling 
computational  frame  with  respect  to  -Z^  plane 
(angle  of  attack  plane) 

Dummy  aerodynamic  roll  angle;  reduced  aerodynamic 
roll  angle;  cp  modulo  where  ^ is  the  repeat  roll 
angle  of  the  aerodynamic  symmetry 

Earth  latitude 
Inertial  latitude 
Target  latitude 

Angular  velocity  of  a mass  particle  within  the 
geometric  boundary  of  the  vehicle  with  respect  to 
rolling  body  frame 

Rotational  speed  of  earth 

dt 

d' 

dt^ 

Transpose  of  matrix 
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I 


Inverse  of  matrix 


Initial  values  unless  otherwise  specified 
Instantaneous  value 

V'alue  of  a variable  when  calculations  are  auto- 
matically terminated 


B.  Fortran  Array  Designations  and  Definitions 


KY  (1) 


KY  (2) 


Run  number 


KY  (3) 


Mode  of  integration-  1 = Runge-Kutta 

2 = Adams-Moulton 

3 = Runge-Kutta  with  truncation 

error  estimate 


KY  (4) 


Number  of  differential  equations;  always  input  19 


KY  (5) 


Always  zero 


KY  (6) 


Always  one 


KY  (7) 


For  no  wind  input  zero;  for  wind  input  one 


KY  (8) 


KY  (9) 


Zero  = no  [C^]  correction;  1 = correct  [E  ] 
Print  frequency  for  output 


KY  (10) 


KY  (11) 


Print  frequency  for  orthogonality  of  [8  ] 

D 

Rotational  symmetry  of  body 


1 = 90°  symmetry 

2 = 180°  symmetry  or  body  of 

revolution 

3 = 360°  symmetry 

4 = 60°  symmetry 


5. 


I 


KY  (12) 


Input  one  to  omit  printing  of  aerodynamic  coefficient 
tables 


KY  (13) 

KY  U4) 

KY  (15) 

KY  (16) 

KY  (17) -KY  (20) 


Y (1)-Y  (9) 

Y (10)-Y(18) 

Y (19)-Y  (27) 

Y (28)-Y  (36) 

Y (37)-Y  (45) 

Y (46) 

I 

I Y (47) 

I 

I 

- Y (48) 

Y (49) 

Y (52) 

t 

I 

t 


Input  one  to  omit  printing  output  formats  4 and  5 

Do  not  imput;  used  in  MAIN  to  call  SETUP  in  PANOL  3 

Input  a number  greater  than  zero  if  the  target 
position  is  input  as  tables  in  feet;  input  zero  if 
a single  value  in  feet;  input  a negative  value  if 
tables  of  latitude  and  longitude  in  radians 

Plot  frequency  of  seeker  field  of  view 

Not  presently  used 


Rotation  matrix,  about  X,  by  d (Subroutine  SETUP) 

Rotation  matrix,  about  Y,  by  a (Subroutine  SETUP) 

Rotation  matrix,  about  X,  by  0-90°  (Subroutine  SETUP) 

Rotation  matrix,  about  Z,  by  y.  (Subroutine  SETUP) 

Product  of  above  matrices  in  subroutine  SETUP; 

[£^]  matrix  in  subroutine  OUT 

Cos(\)/g)  or  cos(vlZj^)  which  are  initially  equal 
(Subroutine  SETUP) 

Sin(v|/_)  or  sin(v|/  ) which  are  initially  equal 
E K 

(Subroutine  SETUP) 


Cos(A^)  (Subroutine  SETUP) 

Sin  (A  ) (Subroutine  SETUP) 

Cos  (vl»  ) or  cos(t  ) , which  are  initially  equal 
£ H 

(Subroutine  SETUP) 
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r*  'S,  i 'r  rWt  Km  . 

- 


I 


i; 

J 


Y (53) 

Sin(t_)  or  sin(x  ) , which  are  initially  equal 

£ r\ 

(Subroutine  SETUP) 

Y (54) 

Cos(y^)  (Subroutine  SETUP) 

Y (55) 

Sin  (y^)  (Subroutine  SETUP) 

Y (56) 

Cos  (a)  (Subroutine  SETUP) 

Y (57) 

Sin  (a)  (Subroutine  SETUP) 

Y (58) 

Cos(0)  (Subroutine  SETUP) 

Y (59) 

Sin(0)  (Subroutine  SETUP) 

Y (66)) 

LE 

Y (67) 

''y 

LE 

Y (68) 

^LE 

Y (69) 

^E 

Y (70) 

^E 

Y (71) 

^E 

Y (72) -Y 

(80) 

[£j^]  array 

Y (81.  . 

.83) 

Cross  velocity  in  computational  frame 

Y (84.  . 

.86) 

Cross  velocity  in  local  frame 

Y (87.  . 

,89) 

Sun  vector  in  local  frame  (PANOL  4) 

Y (90) 

Vv 

BA 
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Y (91) 


Y (92) 


Y (93) 


Y (94) 


Y (95) 


Y (96)-Y  (104) 


Intermediate  storage 


Y (105) 


Y (106) 


Q,  dynamic  pressure 


Y (107)-Y  (109) 


Intermediate  storage 


Y (110) 


Mach  number 


Y (lll)-Y  (116) 


Intermediate  storage 


Y (117) 


Y (118) 


Y (119) 


Y (120)-Y  (130) 


Intermediate  storage 


Y (131) 


Y (132) 


Y (133) 


Y (134) 

Y (135) 

Y (136) 

Y (137) 

Y (138) 


Y (139) 

Y (140) 

Y (141) 

Y (142) 

Y (143) 

Y (144) 

Y (145) 

Y (146) 

Y (147) 

Y (148) 

Y (149)-Y  (162) 

Y (163)-Y  (189) 

Y (190) -Y  (199) 

Y (200) 

Y (201) 


m 

q 

C (originally  yaw  plane  value;  not  now  used) 

"'q 

I 

Intermediate  storage 

F : force  in  the  x,_  direction  (pounds) 

F^;  total  force  in  the  y^^  direction  (pounds) 

F ; total  force  in  the  direction  (pounds) 
z AB 

M ; total  moment  in  the  X direction  (pound-feet) 

M : total  moment  in  the  y,„  direction  (pound-feet) 
y AB 

M ; total  moment  in  the  z . „ direction  (pound-feet) 
z AB 

Intermediate  storage 

Used  in  Subroutine  SETUP  only 

Intermediate  storage 

h^:  terminal  altitude  (feet) 

t^:  terminal  time  (seconds) 
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Y (202) 

Y (203) 

Y (204) 

Y (205) 

Y (206) 

Y (207) 

Y (208) 

Y (209) 

Y (210) 

Y (211) 

Y (212) 

Y (213) 

Y (214) 

Y (215) 

Y (216) 

Y (217) 

Y (218) 


a^;  terminal  angle  of  attack  (degrees) 

A„;  altitude  above  which  atmosphere  is  assumed 
n 

negligible  (feet) 

w-  rotation  rate  of  the  Earth  (radians/ second) 

E 

Error:  zero  to  limit  relative  truncation  error,  minus 
one  to  limit  absolute  truncation  error 

ENE:  fixed  step  if  zero,  truncation  error  level  if 
positive 

At^:  initial  integration  step  size  (seconds) 

t : initial  time  (seconds) 

0 

k;  if  zero,  computational  frame  non-spinning  if  one, 
computational  frame  if  body-fixed 

tii  time  of  main  motor  ignition  (seconds) 

t2:  time  of  main  motor  burnout  (seconds) 

m^:  initial  mass  of  powered  vehicle  (slugs) 

: initial  axial  moment  of  inertia  of  powered  vehicle 
(slug-ft*) 

1 : initial  transverse  moment  of  inertia  of  powered 
^o 

vehicle  (slug-ft*) 
h;  initial  altitude  (feet) 

V^:  initial  speed  with  respect  to  the  Earth  (feet/ sec) 

A^:  azimuth  of  launch  velocity  (degrees) 

Y„:  elevation  angle  of  launch  velocity  (degrees) 
is 
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1 


Y (219) 


0;  direction  with  respect  to  the  Earth  of  plane  of 
angle  of  attack,  see  Figure  5 (degrees) 


Y (220) 

Y (221) 

Y (222) 

Y (223) 

Y (224) 

Y (225) 

Y (226) 

Y (227) 

Y (228) 

Y (229) 

Y (230) 

Y (231) 

Y (232) 

Y (233) 

Y (234) 

Y (235) 

Y (236) 

Y (237) 


f.  initial  longitude  if  geographic  output  required 
(degrees) 

v|/;  initial  latitude  if  geographic  output  required 
(degrees) 

a;  initial  total  angle  of  attack  (degrees) 

0;  initial  direction  of  plane  of  angle  of  attack  in 
computational  frame  (degrees) 

p;  initial  spin  rate  (radians/second) 

q:  initial  pitch  rate  (radians/second) 

r;  initial  yaw  rate  (radians/second) 

V,„:  wind  speed  (feet/ second) 
w 

A^:  wind  direction  azimuth  (degrees) 

R,,:  radius  of  the  Earth  (feet) 

1^:  axial  moment  of  inertia  (slug-feet-squared) 

ly!  transverse  moment  of  inertia  (slug-feet-squared) 

I : transverse  moment  of  inertia  (slug-feet-squared) 
z 

m;  mass  of  vehicle  at  burnout  (slugs) 

A;  vehicle  reference  area  (square  feet) 
d;  vehicle  reference  length  (feet) 

6;  effective  configurational  asymmetry  angle  (radians) 

6 : pitch  asymmetry  angle  (rad  or  deg  depending 
on  c ) 
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Y (238) 

y (239) 

Y (240) 

Y (241) 

Y (242) 

Y (243) 

Y (244) 

Y (245) 

Y (246) 

Y (247) 

Y (248) 

Y (249) 

Y (250) 

Y (251) 


6_;  yaw  asymmetry  angle  (rad  or  deg  depending 

D 

on  C ) 

X ; initial  axial  location  of  entire  vehicle  center 
c 

o 

of  mass  (feet) 

Ci;  thrust  malalignment  in  yaw  plane  (radians) 

Ej;  thrust  malalignment  in  pitch  plane  (radians) 

x^:  axial  location  of  nozzle  exit  in  body  for  thrust 
malalignment  (feet) 

y^;  y-lateral  location  of  nozzle  for  thrust  malalign- 
ment (feet) 

z^;  z-lateral  location  of  nozzle  for  thrust  malalign- 
ment (feet) 

• 

m;  mass  loss  rate  from  propellant  burning 
(slugs/ second) 


m ; initial  mass  of  propellant  (slugs) 


p^:  pressure  of  propellant  gas  at  exit  of  nozzle 
(pounds/ft*) 


A ; nozzle  exit  area  (square  feet) 

6 

V ; velocity  of  propellant  gas  at  exit  (feet/ second) 

0 

Xi^;  initial  position  of  front  face  of  propellant  grain 
(feet) 


X2  : initial  position  of  center  of  mass  of  propellant 
0 

grain  (feet) 


Y (252) 

Y (253) 

Y (254) 

Y (255) 

Y (256) 

Y (257) 

Y (258) 

Y (259) 

Y (260) 

Y (261) 

Y (262) 

Y (263) 

Y (366) 

Y (367) 

Y (368) 


Xj  : initial  position  of  rear  face  of  propellant  grain 
o 

(feet) 

X4;  longitudinal  position  of  nozzle  exit  plane  for 
jet  damping  (feet) 

rj;  outer  radius  of  propellant  grain  (feet) 

T2  ; initial  inner  radius  of  propellant  grain  (feet) 
o 

p^;  density  of  propellant  grain  (slugs/cubic  feet) 

x^',  axial  shift  of  burnout  center  of  mass  due  to 
EC 

vehicle  unbalance  (inches) 

y : lateral  shift  of  burnout  center  of  mass  due  to 
vehicle  unbalance  (inches) 

k;  azimuthal  orientation  of  lateral  center  of  mass 
position  due  to  vehicle  unbalance  (degrees) 

Specific  impulse  (lb  sec/lb) 

Hj:  first  Euler  angle  in  transformation  from  com- 
putational frame  to  principal  axis  frame , about 
z axis  degrees) 

H2;  second  Euler  angle,  about  y'  axis  (degrees) 
Hj;  third  Euler  angle,  about  x"  axis  (degrees) 


Aerodynamic  torques  due  to  eccentric 
center  of  mass  position 
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Y (369) 

Y (370) 

Y (371) 

Y (373) 

Y (373) 

Y (374) 


Y (375)-Y  (383) 

Y (384)-Y  (392) 

Y (393)-Y  (401) 

Y (402) -Y  (410) 

Y Y (413) 

(416) 

Y (417)-Y  (425) 

Y (426)-Y  (434) 

Y (435)-Y  (437) 

Y (438) 

Y (439) 

Y (440) -Y  (448) 


Y (449)-Y  (451) 


Accelerations  in  computational  frame  for 
output  tape 


Velocity  components  in  computational  frame  for 
output  tape 

Available  for  use 

Inertia  tensor  for  principal  axes 

Intermediate  storage 

Inertia  tensor  for  skewed  axes 

Angular  momentum  H 

Omega  time  H 

Transpose  of  inertia  tensor 
Inverse  of  inertia  tensor 
Final  sum  of  torques 

Angle  between  sun  vector  and  body  axis  (SIGMA) 
THETA 

Tensor  RTEN  to  translate  inertia  tensor  to  off-axis 
center  of  mass 

Missile  axis  components  in  Local  Frame;  Panol  4 only 
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Y (452)-Y  (455) 

Y (456) 

Y (457) 

Y (458) 

Y (460)-Y  (468) 

Y (470)-Y  (478) 

Y (480)-Y  (488) 

Y (490) -Y  (498) 

I 

Y (601) -Y  (3707) 


q,  r.  q,  r in  body-fixed  (rotating  frame 
Spin  number  (non-dimensional  spin) 
Beta  in  non-spinning  triad  for  plotting 
Alpha  in  non-spinning  triad  for  plotting 

Inertia  tensors  in  PANOL  3 only 

Table  storage  for  output 
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STATE  VARIABLE-C  ARRAY 


c (1) 

m 

C (2) 

• 

C (3) 

^R 

C (4) 

^R 

C (5) 

^R 

C (6) 

"r 

C (7) 

P 

C(8) 

q 

C(9) 

r 

C (10) 

C (11) 

C (12) 

C (13) 

[«b1m 

C (14) 

C (15) 

[fifils* 

C (16) 

r 


c 

c 

c 

c 

c 


(17) 

(18) 

3 3 

(19) 

cp' 

(20)-C  (21) 

Not  presently  used 

(22) 

Error  control 

C (23) 


Not  presently  used 


STATE  VARIABLE  DERIVATIVE;  D ARRAY 


D U) 

D(2) 

• m 

D (3) 

^R 

D (4) 

^R 

D (5) 

• 

^R 

D (6) 

• 

"r 

D (7) 

• 

P 

D (8) 

• 

q 

D (9) 

% 

r 

D (10) 

D (11) 

D (12) 

D (13) 

t^Bl 

D (14) 

t«B^ 

D (15) 

t^B^ 

D (16) 

f^*Bl 

• 

D (17) 

• 

D (18) 

• 

D (19) 

9' 
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